A Threestepped Coordinated Level Set Segmentation Method for
Identifying Atherosclerotic Plaques on MR-Images
Oliver Gloger*, Matthias Ehrhardt**, Thore Dietrich***, Olaf Hellwich, Kristof Graf***,
Eike Nagel****
Berlin University of Technology (TUB), Computer Vision and Remote Sensing
Franklinstraße 28/29, Sekretariat FR3-1, 10587 Berlin, Germany
Tel.: ++49-30-314-73109 , Fax: ++49-30-314-21104,
berlin.de, [email protected], [email protected]
Abstract
In this work we propose an adapted level set segmentation technique for the recognition of
atherosclerotic plaque tissue on magnetic resonance images. The images are 2dimensional
cross-sectional images and show different profiles from ex-vivo human vessels with high
variability in vessel shape. We used a curvature based anisotropic diffusion technique to
denoise the magnetic resonance images.
The segmentation technique is subdivided into three level set steps. Hence, the result of every
phase serves as constructive knowledge for the next level set step. By analyzing and
combining carefully all available channel information during the first and second step we are
capable to delineate exactly the vessel walls by using and adapting two well-known level set
segmentation techniques.
The third step controls an enclosing level set which separates the plaque patterns from healthy
media tissue. In this step we introduce a local weighting concept to consider intensity
information for conspicuous plaque patterns. Furthermore, we propose the introduction of a
maximal shrinking distance for the third level set in the vessel wall and compare the results of
the local weighting algorithm with and without the concept of the maximal shrinking distance.
The incorporation of locally weighted intensity information into the level set method allows
the algorithm to automatically distinguish plaque from healthy media tissue. The knowledge
of the maximal shrinking distance can improve the segmentation results and enables to
delineate tissue areas where plaque is most likely.
Keywords: level set segmentation, active contours, medical image segmentation, anisotropic
diffusion, atherosclerotic plaques, Canny edges
* Chair for Computer Vision and Remote Sensing (TU Berlin) /
Ernst Moritz Arndt University Greifswald, Institute for Community Medicine, Walter-Rathenau-Str. 48,
17475 Greifswald
** Weierstrass Institute for Applied Analysis and Stochastics, Mohrenstr.39, 10117 Berlin, Germany
*** German Heart Center Berlin, Augustenburger Platz 1, 13353 Berlin
**** King's College London, St. Thomas’ Hospital, London SE1 7EH
1
1. Introduction and Motivation
Atherosclerosis is a serious disease affecting the arterial blood vessels. It is a continuous
process in which certain substances (e.g. Cholesterol) accumulate within the vessel wall. This
accumulation may cause a slow reduction of the vessel lumen causing a deterioration of blood
flow during exertion, or – in the case of a “vulnerable” plaque - rupture suddenly causing a
closure of the vessel, which leads to a myocardial infarction (heart attack) or stroke. Thus, it is
of most importance to discriminate between stable and unstable plaque and predict potential
plaque rupture.
With magnetic resonance (MR) imaging (MRI) it is possible to visualize all parts of the body
– including the heart and the vessels - with high spatial resolution (typical values for in-plane
resolutions in clinical research are about 400 µm x 400 µm for human plaques in coronary
arteries and about 500 µm x 500 µm for aorta plaques). It is possible to visualize plaque in
larger vessels, such as the aorta and the carotid arteries with state of the art MR scanners and
to characterize plaque components. Methods for automatic plaque segmentation of carotid
plaque have been proposed [1] and allow for serial assessment of plaque size and its lipid
core. Similarly vessel wall imaging of the coronary arteries is possible and first clinical data is
available [2, 3, 4]. However, the rapid motion of the coronary arteries in combination with the
small vessel size has made plaque characterization in the coronary arteries impossible.
Potential future developments will improve spatial resolution, overcome motion artifacts by
motion compensation algorithms [5, 6] or specifically enhance those components of a plaque
which cause plaque rupture using specific contrast agents [7, 8].
Due to prospective improvements of MR-image resolution we aim to test modern
segmentation techniques for higher in-plane resolutions of 79 µm x 79 µm, which is of course
five to six times higher than usual. Informations about segmentation success for this type of
MR-images should help to estimate the benefit of level set segmentation techniques for plaque
images in general. If we expect higher resolutions for aorta or coronary arteries in future, this
contribution represents a pioneer work for automatic segmentation of prospective image
resolutions.
Primarily, we applied pixel-oriented classification algorithms (i.e. k-means) and low-level
segmentation approaches (i.e. watershed-segmentation, region growing, presegmentation with
scale space concepts [9]) for plaque identification, but stated that suchlike approaches failed
in most circumstances. To facilitate automatic image analysis and to reduce observer
variability we propose a new algorithm for the assessment of plaque images. In contrast to
previous methods with probabilistic approaches we provide a technique that subdivides
healthy media tissue from plaque with a closed contour into two resulting subparts. This
bordering closed contour guarantees an important local closeness to distinguish between
plaque and healthy media tissue during the segmentation process.
Atherosclerotic plaques cannot be differentiated on MR-images by the information of gray
values alone, but have an associated regional context. We know however, that atherosclerotic
plaques are located within or adjacent to the arterial walls. Consequently, we are interested in
a segmentation technique that – in a first step- can segment the arterial vessels and then – in a
second step- can recognize relevant plaque structures in the vessel walls. Since vessel
boundaries in 2D cross-sections must present closed loops, we used a segmentation technique
based on active contours as an appropriate starting point.
2
The model of active contours has been introduced by Kass et. al. [10] and represents an
explicit form for describing the motion of a closed contour. The algorithm is based on the
principle, that a contour continues to move towards the boundaries of an object. This process
can be represented by a minimization of an explicit energy function which includes two
terms: one for controlling the smoothness of the contour and one to influence the attraction of
the contour to the object boundaries. This parametric model has been frequently used and
extended for segmentation problems in (bio)medical image segmentation [11, 12, 13, 14, 15,
16].
The classical model of snakes and active contours yields a significant disadvantage of
managing certain topological changes during the contour evolution process. Operations like
splitting and merging are not possible to represent with this classical approach. Unfortunately,
those flexible curve propagations could be of great importance during the segmentation
process of detailed medical shapes especially for the delineation of the outer media border and
the recognition of plaque (which we show in detail in chapter 4.2). Caselles et al. [17]
associated the classical active contour model with a geodesic flow model which is represented
by a partial differential equation based on a mean curvature flow. The energy minimizing
problem can be solved by finding the geodesic contour in Riemann space derived from the
image. Solving this model with a level set approach provides the contour with the ability of
managing the desired topological changes. The contour is now represented implicitly through
the level set equation and yields additionally the advantage of having fewer parameters to
guide the propagation of the contour. Parameters used in the level set propagation are
geometric measures which are independent of the parametrization of the curve. So they are
also called geometric active contours where the classical snakes refer to parametric active
contours.
In combination with variational techniques the level set method has often been applied for
finding the minima of energy functionals, which aim to solve special segmentation problems
for modern computer vision applications. In [18] a dynamical framework is proposed that
uses variational techniques derived from a Bayesian model and is able to segment scalar and
vector-valued images. In contrast to the edge based energy functionals [17], region based
energy functionals [19, 20] are independent of edge information making their application for
medical images very promising. However, variational formulations that use either edge based
terms or region based terms have the disadvantage that their derived level set propagation
stops in local minima. Recently, some approaches try to overcome those problems. In [21] a
local term is incorporated in a region based energy functional [19] that minimizes the local
contrast variances on a narrow band inside and outside of the contour. Lankton et al. [22]
combine the advantages of region based [19] and edge based energy functionals [17] by
deriving a geodesic energy from local regions around the contour which results in a hybrid
curve evolution.
The exclusive application of region based and edge based energy functionals is not sufficient
for atherosclerotic MR-images due to their vessel complexity that show high shape and
intensity variances and lacks of local edge information. In this work we propose a level set
technique that uses relevant informations during propagation to overcome those challenges
and is capable for plaque delineation on a variety of MR-images.
In the following Chapters we explain this new approach for plaque segmentation successively.
In Chapter 2 we describe the image data we work with and explain how they have been
created. In Chapter 3 we describe the denoising method we used for the MR-images as
preprocessing step in our algorithm. In Chapter 4.1 we give a short summary of the theoretical
3
aspects of the level set method which are crucial for understanding the subsequential
concepts. Chapter 4.2 describes our threestepped level set method in detail and shows the
results of the segmentation process of the inner and outer media borders. Furthermore, this
Chapter introduces the proposed concepts for the local weighting of conspicuous vessel wall
areas and the healthy media border. Chapter 5 illustrates and compares the results for the
proposed concepts. Finally, in Chapter 6 we give a summary and a conclusion of this work
and summarize our ideas for future work.
2. MR Data Acquisition
Human arterial specimens from the groin (iliac region) were obtained during vessel
replacement surgery. After freezing at -80°C, each specimen was warmed up to 30°C in a
phosphate buffered solution. Then it was placed in a 15 ml tube which was then filled with
1% agarose solution. The specimens were measured in a 7 Tesla MR-Tomograph (Bruker
Germany, PharmaScan). A whole body mouse coil (diameter of 38 mm) was utilized to obtain
a 79 µm x 79 µm in-plane resolution and a 110 µm slice thickness. Scan time for four
contrasts (PD, T1, T1fatsat, T2) was about 16 hours. The method is described in more detail
in [23]. Figure 1 shows an imaging example of 4 channels each with a 256 x 256 matrix size.
With those different channels we have the possibility to emphasize special tissue patterns in a
certain range.
PD
T1
muscular vessel
wall (media)
chronic plaque
young plaque
(neointima
formation) Lumen
T1fatsat
T2
calcification
Fig.1: MR-images of arterial profiles quantified with 4 different types of measurement (proton-density, T1, T1-
fatsaturated, T2). The important parts of the arterial profiles are inscribed to the images above.
4
3. Denoising of MR-images
Although certain relevant plaque patterns are well visible the images exhibit regions with a
low signal-to-noise ratio (SNR). The level set propagation is influenced by high image
gradients. To ensure a suitable level set propagation we have to denoise those images
effectively without loosing too much detail information. Blurring the images with a Gaussian
kernel causes an additional dislocation and blurring effect for important edges. Those edges
are indispensable information for the stopping decisions of the level set propagation.
Consequently, we decided to apply edge-preserving algorithms based on anisotropic
diffusion.
A frequently referenced non-linear diffusion model was introduced by Perona and Malik [24].
In [25, 26] two important extensions of the Perona-Malik model have been proposed. Wei
[25] suggests a generalisation of the Perona-Malik equation by introducing an edge enhancing
functional and using local statistics for more effective edge preservation and noise removal.
Another extension of the Perona-Malik model (with representing the underlying image) has
been developed by Whitaker and Xue [26], which they call a modified curvature diffusion
equation (MCDE):
u
(( ) )
uu
udivg u
tu
∂∇
=∇ ⋅ ∇ ⋅
∂∇
(1)
They propose it as a level-set analogue of the anisotropic diffusion model of Perona and
Malik and also emphasize its superiority compared to the Perona-Malik model. Their curve
evolution is based upon the principles of the self-snakes, which have been introduced by
Sapiro [27].
Due to the fact that we use denoising techniques as a preprocessing step for our work, we
chose the implementation of the ITK development framework [28] to denoise MR images
using the MCDE model. The diffusion process can be controlled by the number of iterations
and a conductance value, which influences the conductance term of the diffusion process. In
the ITK the conductance term [28] is used as a function that can suppress diffusion at higher
image gradients and preserve edges by
()
2
2
exp 2
u
gu k
⎛⎞
∇
∇
=⎜−
⎜⎟
⎝⎠
⎟
(2)
and can be steered by the conductance value k. We used 175 iterations for time steps of 0.125
seconds and denoised the images with conductance values in the range of 1.0 to 2.5.
5
4. Plaque Segmentation with the Level Set Method
4.1. Description of the Level Set Method
The level set method traces back to the work of Dervieux and Thomasset [29] and has been
introduced for the purpose of curve evolution by Osher and Sethian in 1988 [30]. Level set
methods gained increasing interest in the field of fluid dynamics and material sciences and
became also popular in computer vision fields mainly in medical image processing. (A couple
of level set applications are described in [31, 32, 33]). The main idea in the level set concept
is describing the subsequent motion of an implicitly given closed contour via a time-
dependent partial differential equation (PDE). Instead of dislocating the contour explicitly,
what has to be arranged for parametric active contours, we embed the initial location of the
two-dimensional front as a zero level set of a three-dimensional function (, ,)
x
yt
φ
. This
function is a signed distance function, which provides for every location in the image its
Euclidean distance to the zero level set. Using this definition the resulting movement of the
curve is directly related to the temporal dislocation of the zero level set (Fig.2).
Fig.2: Embedding a closed contour as a zero level set into a 3D time-dependent function. The 3 resulting zero
level sets are shown in 2D with their corresponding values.
If we calculate the time-dependent alteration of the signed distance function for certain time
steps, we can follow the time-dependent displacement of the curve as the zero level set of the
signed distance function. The PDE can be solved by using explicit or implicit methods for
time discretization. By using explicit methods we have to ensure numerical stability, which
can be realized by applying relatively small time steps for the time discretization that fulfill
the condition of Courant-Friedrichs-Lewy (CFL condition).
Segmentation techniques using level set propagation and considering image properties (i.e.
edge information) and curve properties (i.e. curvature) were suggested by Malladi et. al. [34]
and Caselles et. al. [17]. Stopping terms for implicit active contours with mean curvature
motion have been introduced in [35]. The authors in [17, 34, 35] propose stopping terms,
which are based on image gradients and are represented by the factors:
0
1
() , 1
1(,)(,)
p
o
gu p
Gxyuxy
σ
=≥
+∇ ∗ (3),
0
0
1
() exp( ( , ) )
gu Guxy
σ
=∇∗ (4)
6
with 0
u as the original image and G
σ
as the Gaussian function with standard deviation
σ
.
The curvature can be expressed by the divergence of the normal, which can be simplified to
the Laplacian, if we ensure the condition 1
φ
∇
=. This can be guaranteed by a reinitialization
process which aims to renew the signed distances after every time step.
We apply a segmentation technique for the outer media border, which relies upon the
principle of Malladi et. al. [34], who introduced curvature and image gradients into the level
set propagation via
(
)
0
() 0
tAG
gu F F
φφ
+⋅+⋅∇=
(5)
where the force represents the force in normal direction for the level set propagation.
Depending on the sign of the distance values inside and outside the curve, the sign of can
be used to shrink or expand the evolving curve. The other term stands for the geometry of
the front such like the local curvature. With the model (5) we can influence a zero level set to
stop at the desired boundaries of the object and penalize its local curvature. The penalization
of the curvature is important for medical objects, which in general exhibit smoothed
boundaries. However, this model requires edge information for the whole boundary of the
object which is to be segmented. For medical images the problem of missing edges is well-
known, so the application of (5) can lead to the leakage of the zero level set through the gaps
caused by missing edges. A possibility to resolve this problem is the geodesic active contour
model proposed by Caselles et al. [17]. They show that finding the geodesic curve in
Riemannian space is equivalent to solving the classical energy model of the active contours.
The geodesic active contour model (6) yields an additional attraction term with respect to (5):
A
F
A
F
G
F
00
()( ) ()
tAG
gu F F gu
φ
φφ
=⋅+⋅∇+∇⋅∇
(6)
The scalar product of the attraction term g
φ
∇
⋅∇ attracts the evolving curve to the boundary,
if the curve is near the boundary. The term can even prevent the leakage of the evolving curve
through weak edges and small gaps.
4.2 The Threestepped Level Set Segmentation
The goal of our method is the segmentation of plaque from MR images measured in [23]. The
media can exhibit a variety of plaque structures, however those formations arise adjacent or
inside the vessel wall. We aim to segment the media together with its plaque tissue and after
that delineating healthy tissue from plaque tissue inside the vessel wall.
For the segmentation we used the level set method. The process of our threestepped plaque
segmentation can be subdivided in three phases. Every phase uses a single level set approach
and the results of every phase serve as additional knowledge for the next level set step.
Primarily, we apply an inner zero level set to separate the lumen from the inner media border.
After that we evolve the final contour of the first level set phase until it reaches the outer
boundary of the media.
7
The first two steps ensure the segmentation of the media. The most challenging task is the
differentiation of plaque from healthy media tissue inside the extracted media. Fortunately,
we know that atherosclerosis is a continuous process beginning adjacent to the inner media
border. Neointimal formations accumulate at the inner media border and chronic plaque
develops successively inside the media. Subsequently, we exploit the closed contour of the
zero level set to delineate the healthy media tissue from plaque with a shrinking zero level set
starting from the outer media border. This coordinated process will be explained in the
following three sections in more detail.
4.2.1. Detecting the Inner Border of the Media
We expect the outer border to be approximately circle-shaped. In our images of surgical
specimen this is not fulfilled due to deformations during surgery and further preparation. The
inner border of the vessel can have more variability, but will always form a closed shape with
rounded inner edges. In our considerations we exclude vessels with a completely closed
lumen. Consequently, the initial seed point was placed in the center of the arterial profile for
every image.
On T1 and T1fatsat images exist more contrast between the media and its adjacent regions in
comparison with the corresponding PD and T2 images. To preserve the contrast for the media
boundaries we denoised the T1 and T1fatsat images by using a lower conductance value of
1.0. Due to possible artifacts and signal descents during measurement the MR images can
loose information and miss edge magnitude. Executing the algorithm for only one channel
image could cause leakage effects for the zero level set. We aim to compensate this effect by
averaging the intensities from two measurements. Although contrast can be deteriorated in
places, we collect sufficient edge magnitude to segment the media by averaging the denoised
T1 and T1fatSat images.
For guiding the first level set to the inner media border we used the geodesic active contour
model [4] and used (6) in combination with (3) (p=2). Starting from an initial point inside the
lumen, the first zero level set terminates its propagation exactly at the inner media borders.
Near the inner media border the attraction term attracts the curve to the exact edges of the
border and furthermore avoids the leakage into the vessel wall region (s. Fig.3). On its way to
the inner media border the zero level set encircles those edges, that result from
inhomogeneities inside the lumen.
8
Fig.3: Curve evolution of the geodesic active contour model converging to the inner media border. Three time
steps with the zero level set (green curve) are shown . Left: zero level set on edge image relying upon the
stopping term (g). Middle: Attraction force ( g
φ
∇
⋅∇ ) from the geodesic active contour model. Right:
propagating curve on denoised image. (We used propagation speed 6F
=
and curvature speed ). 2F=
4.2.2. Detecting the Outer Border of the Media
It is not possible to segment the outer media border with a shrinking level set starting beyond
the media. There could exist parts from other tissues outside of the vessel wall which could
unintentionally stop a shrinking zero level set before it reaches the outer media border. After
delineating the inner media border there is not reliable knowledge about the extent of plaque
and media tissues. Subsequently, it is precarious to estimate an initial position for a shrinking
level set beyond the media based upon the location from the delineated inner media border.
Hence, we utilize the position of the inner media border as determined in the first step as
starting zero level set for the next segmentation step.
Due to higher inhomogeneities of plaque tissue inside the vessel wall there exist stronger
edges than in the lumen region. For the further propagation we stated that the geodesic active
contour model can be attracted by strong edges inside the vessel wall and stop the propagation
unintentionally. Hence, we extracted the relevant edges with the Canny algorithm [36] and
used Canny edges as stopping criterion. (For the hysteresis in the Canny algorithm we chose
0.05 as higher threshold and 0.02 as lower threshold for all images). A similar approach to
combine the level set propagation with the Canny algorithm has been proposed by Quin [37],
who use the maximum along the gradient directions to determine the edge points. For the
second level set step we applied the concept of [34] and used (5) to guide the second level set
to the outer media border. For using the flexible motion character of the level sets we set the
curvature speed equal to zero ( 0
G
F=), to allow the zero level set to circulate around Canny
edges inside the vessel wall. Instead of using (3) or (4) we used Canny edges as stopping term
9
0, (,) _
(, ) 1,
if x y is Canny edge
canny x y else
⎧
=⎨
⎩ (7)
To overcome the barricades of the edge of the inner media border we let the zero level set
pass by using one iteration without applying the stopping term. Similarly to the first step, we
apply the second segmentation step to the average image of the denoised T1 and T1fatsat
images. To exploit the edge enhancement of the anisotropic diffusion we used conductance
values of 2.5 to denoise the images. The denoising method preserves all Canny edges of the
outer media border and attenuates edges inside the vessel wall.
Fig. 4: Approaching the second level set (blue curve) to the outer boundary of the arterial vessel. The selected 4
time steps are listed pairwise. The final zero level set of the inner vessel border is shown in green colour for
every time step in the right image. The circumfluent propagation (propagation speeds of 2 ) along Canny
edges can be seen in the left image of every time step.
A
F=
After this step the whole media is delineated which provides a comfortable initial point for
further plaque recognition, because atherosclerotic plaque can only be located adjacent to the
inner media border or in the media itself.
4.2.3. Contouring of Atherosclerotic Plaque Patterns in the Media
In contrast to the first and second level set step we need as much information from plaque
tissue as possible for the third level set step. Plaque tissue is enhanced in PD and T2 channel
images. However, in [23] it was argued that T2 weighted images exhibit most plaque contrast,
but also in PD weighted images plaque is enhanced. For collecting as much plaque
information as possible we denoised the PD and T2 images with a conductance term of 1.0
and used the average image for the third level set step.
To present our results we provide the numbers for the 12 different specimens which were used
in [23] and show their T2-measurements in Fig. 5.
10
Fig. 5: T2 images of 8 samples with recognizable plaque (image numbers: 9, 11-15, 19, 20) and 4 samples
without plaque tissue (image numbers: 10, 16-18).
The plaque images do not show suitable edge information with respect to the relevant plaque
tissue. High gradient magnitudes inside the vessel wall influence the zero level set to stop too
early, whereas Canny edges exhibit gaps. Consequently, we need a further concept, that
supports plaque segmentation inside the vessel wall.
4.2.3.1. Penalizing Grayvalue Deviances using a Local Weighting Concept
Regarding the plaque images from Fig. 5 it gets obvious that plaque patterns attract attention
also with their distribution of intensity information. Neointimal formations yield higher
intensities whereas chronic plaque is darker in comparison with the adjacent healthy media
tissue. Thus, we decided to penalize striking intensities additionally to Canny edges during
the level set propagation.
To obtain a measure for the intensity distribution of healthy media tissue, we consider the
intensities in a small stripe inside the media adjacent to the outer media border. We disregard
the case of a media where chronic plaque is tangent to the outer media border, because all our
analyzed images show at least such a small stripe without inflamed media tissue. We apply a
11
shrinking level set with negative propagation speed (-3) for two iterations to delineate a small
stripe between the actual zero level set and the outer media border. However, it is precarious
to consider the intensity distribution of the whole stripe to distinguish plaque and healthy
media. During MR measurement can arise artifacts and descents of signal intensity in an
arbitrary direction. Due to those technical features, it is not predictable, if an intensity at a
certain location on the media does represent the same tissue property of another media
location. Subsequently, the intensities in the whole stripe are not reliable measure for healthy
media intensities. To overcome this challenge we have to incorporate locally weighted
intensity information from the stripe into the level set equation.
For every pixel on the outer media border (we call it the reference point) we calculate a
weighted average value and an accordant weighted standard deviation. All pixels inside the
stripe contribute to the weighted average value and weighted standard deviation of a reference
point with their reference distance. The reference distance from a certain point inside the
stripe is the length on the outer border between its closest point on the outer border and the
reference point. Due to the curvature values for the outer media border, the reference distance
ensures an appropriate distance measure for our weighting method. For the construction of the
weighted average and the weighted standard deviation for every reference point we provide
the weights of the points inside the stripe with a Gaussian weighting function:
2
2
1(
( , ) exp( )
2
2
dxy
wxy
σ
σπ
=⋅−
,)
(8)
where drepresents the explained reference distance,
σ
the selectable standard deviation and
(, )
x
ya position of a point on the outer stripe. The weighting procedure can be demonstrated
with Fig. 6.
Fig. 6: Local weighting relative to a certain reference point and a selected standard deviation. The weighted
values for the pixels inside the stripe based on a selected standard deviation are shown with a thin Gaussian-
shaped wall. The reference point and the points on its perpendicular in this stripe have the highest weighted
values for the local weighting process.
Choosing a very small standard deviation results in weighting the intensities inside the stripe
more locally. By selecting a high standard deviation for the Gaussian weighting function we
will integrate more information from the stripe to make a local decision for conspicuous
regions. Increasing the standard deviation will converge to a uniformly distributed weighting
for all points inside the stripe.
With the local weighting method we determine the dispersion range for healthy media tissue.
12
The weighted average value (, )mxy for every reference point on the outer media border is
given by
1
1
(, ) (, )
(, )
(, )
N
ii ii
iN
ii
i
wx y Ix y
mxy wx y
=
=
⋅
=∑
∑ (9)
where i should be an index referencing all points i
p
with location (, )
ii
x
y in the outer stripe,
(, )
ii
I
xy the intensity of (, )
ii
x
y and N the total number of points in this stripe. Its
corresponding weighted standard deviation is given by:
2
1
1
(, )((, ) (,))
(, )
(, )
N
ii ii
iN
ii
i
wx y Ix y mxy
std x y wx y
=
=
⋅−
=∑
∑ (10)
Every pixel in the media, which is not located inside the stripe, has a closest point (reference
point) on the outer media border. Hence, we compare every pixel intensity in the media,
whose pixel is not located inside the stripe, with the dispersion range of its reference point. If
the intensity falls inside the range of the weighted standard deviation about the weighted
average of the reference point, then this pixel belongs to healthy media tissue. If this is not the
case then this intensity is locally striking and belongs to plaque tissue. Consequently, we
penalize those intensities exponentially with their distance to the dispersion range of the
weighted standard deviation. We express this intensity penalty with:
1,(,)(
(, )
exp( ( ( , ) ( , ) ( , ))) ,
rr rr
rr rr
,)(,)
I
xy mx y stdx y
hxy Ixy mx y stdx y else
⎧−<
⎪
=⎨−− −
⎪
⎩
(11)
In this context the intensity of an evaluated pixel location (, )
I
xy is compared with the
weighted average (
rr
mxy)
and weighted standard deviation (, )
rr
std x y resp. its reference
point (, )
rr
x
y.
We incorporate the penalty of (11) into the level set equation. Furthermore, we know the pixel
locations of the lumen from the first level set step and avoid a possible leakage into the lumen
area with the stopping criterion of
0, (,)
(, ) 1,
if x y lumen
lumen x y else
∈
⎧
=⎨
⎩ (12)
Subsequently, our level set equation results in
A
canny h lumen F
t
φ
φ
∂=⋅⋅⋅⋅∇
∂ (13)
Due to the fact that edges calculated by Gaussian convolution are very broad, they do not
delineate plaque tissue precisely. Canny edges inside the media exhibit gaps, where the zero
level set can break through. Furthermore, Canny edges can divide media tissue that is not
13
relevant for plaque recognition. With the local weighting concept we introduce a measure that
is independent from edge information and regard local grayvalue distributions which are
salient for a human observer. These salient areas are also important for cardiologists to
distinguish plaque tissue from healthy media tissue. Hence, the local weighting approach can
support Canny edges to delineate plaque tissue. Canny edges which do not delineate salient
areas, are encircled by the zero level set.
It should be mentioned that the segmentation results of the local weighting concept can be
steered by the selection of the standard deviation σ in (8). Depending on the circumstances of
the measurement and the interests of the segmentation, it could be important to incorporate
intensities of the stripe which are very close to the reference point. Our purpose was to be
independent of the signal descent of the intensities during measurement, so we used values for
sigma in the range of 10 to 20. Due to the dependence of the outer stripe for the local
weighting it should be clarified that the segmentation results are also dependent on the results
of the first and second level set step. If the media was not segmented correctly, the outer stripe
could contain intensities of other tissue, which could falsify the local weighting results. With
this concept we take into account the shape of the outer border providing the distant values
along the level sets.
4.2.3.2. Introduction of the Maximal Shrinking Distance in the Media
Finally, we use additional information for supporting the segmentation process of plaque. On
most available MR images, that we could apply for this algorithm, there exists at least one
position at the inner media border which does not show plaque accumulation. The closest
Euclidean distance from such a location to the outer media border serves as an estimate of the
healthy media thickness. All tissue located within the outer media border that was not
identified as lumen and has a higher Euclidean distance to the outer media border than the
healthy media thickness delineates a region where plaque is most likely.
A healthy media without plaque accumulation shows very small alterations of thickness.
Hence, we cannot delineate an exact healthy media thickness due to those small thickness
alterations. On plaque images, however, we can delineate the most likely plaque region which
serves as a good approximation for a healthy media border. We use this healthy media border
as a limitation for the shrinking zero level set inside the media. The distance in normal
direction from the outer media border to this healthy media border is called the maximal
shrinking distance (MSD) for the third zero level set. This distance measure is independent of
the vessel wall thickness and can help to avoid the leakage of the zero level set into plaque
regions, which are not penalized by the local weighting concept or delineated by Canny
edges.
On the other hand, if the whole inner media border is coated with neointimal plaque
formations (e.g. image numbers 13 and 15), the smallest Euclidean distance between the outer
and inner media border serves as a helpful measure to delineate most likely plaque regions,
which are located farer away from the outer media border than this distance. Here this
delineating border cannot be regarded as an estimation for healthy media tissue. But it serves
as a maximal shrinking distance (MSD) for the third zero level set, because media regions
beyond this border must be neointimal plaque formations, if we neglect small vessel wall
alterations.
14
For every location in the arterial profile we calculate its smallest Euclidean distance
(dist_to_outer_border in (14)) to the outer media border and compare it with the MSD. We
formulate this as a stopping criterion for the level set propagation in (14):
0, __ _ (,)
(, ) 1,
if dist to outer border x y MSD
smallerMSD x y else
>
⎧
=⎨
⎩ (14)
Due to the fact, that the maximal shrinking distance depends on the most narrow position of
the vessel wall and orientates itself to the outer media border, it is not necessary to regard the
lumen positions given by (12) anymore. We incorporate the knowledge about the MSD into
the level set equation and exchange (13) by (15)
A
canny h smallerMSD F
t
φ
φ
∂
=
⋅⋅ ⋅ ⋅∇
∂ (15)
5. Results
5.1 Results of the local weighting concept
We applied the third level set step according to (13) and list in Fig. 9 the results for the central
of the 3 slices from the analyzed specimens. To compare the quality of the results we regard
expert knowledge about the underlying plaque tissue. For measuring purposes the specimens
have been fixed with agarose gel, which has been manually removed for further histological
examinations. During the histological preparation the vessel form has been extremely altered,
so that a registration between histological images and MR-images is a very complicated issue
(Fig. 7). Furthermore, histological images and MR-images exhibit different slice thicknesses
(110 µm for MR images but 5-10µm for histological images) and due to partial volume
effects an automatic registration between those types of images is awkward.
11 13
vessel
wall
vessel
wall
plaque
tissue
plaque
tissue
Fig. 7: Histological images for image no. 11 (left) and 13 (right). Plaque tissue (neointimal formation) is shown
adjacent to the inner media border. The vessel shapes have been modified during histological preparations.
Reliable knowledge for plaque occurrence has been introduced by experts after comparing
histological images with corresponding MR images. Expert knowledge about plaque regions
is shown in Fig. 8, 9 and 10 (right column) enclosed by blue curves.
15
Fig. 8: Final result of the level set propagation with locally weighted deviant areas for image numbers 9-14. Left: local
deviant areas with Canny edges and segmented lumen region. Middle: Final zero level set after third level set step (red) and
inner media border from first level set step (green). Right: Comparison of result from third level set step (red) with manually
evaluated expert knowledge of plaque tissue (blue) (Standard deviation of the weighting function in the range of 10-20).
16
Fig. 9: Final result of the level set propagation with locally weighted deviant areas for image numbers 15-20. Left: local
deviant areas with Canny edges and segmented lumen region. Middle: Final zero level set after third level set step (red) and
inner media border from first level set step (green). Right: Comparison of result from third level set step (red) with manually
evaluated expert knowledge of plaque tissue (blue) (Standard deviation of the weighting function in the range of 10-20).
17
If Canny edges are connected with delineated plaque the zero level set will stop there. Hence,
in some cases, the final closed contour obtains convexities with a thickness of one pixel.
Those convexities do not represent plaque. A final mean curvature motion with a curvature
penalty of the level sets cannot separate the contour from those convexities due to the closed
connection. We know, that Canny edges have a thickness of one pixel [36]. We exploit this
knowledge and remove coordinates which occur twice in the final zero level set and get rid of
the convexities successfully. The combination success of Canny edges and local deviant areas
to delineate plaque from healthy media tissue gets obvious. Missing edge information is
compensated by local deviant plaque areas. On the other hand, plaque intensities which were
not evaluated as local deviant, are delineated by Canny edges. Furthermore, Canny edges
which do not delineate plaque are encircled by the zero level set.
The image numbers 17 and 18 do not exhibit plaque tissue (Fig. 9). There are no conspicuous
grayvalue areas inside the vessel wall that could be penalized by our algorithm, hence, for
image no. 18 the algorithm does not find plaque. The small light artifact in image number 17
is recognized by the local weighting concept, but has no connection to the lumen. Hence, the
final zero level set excludes this artifact.
The denoising method removes very thin plaque patterns, hence we loose edges and
informations about local deviant areas to stop the level set propagation. However, we regard
those thin plaque patterns as not relevant for our segmentation results. The very small stripe
of plaque tissue at the inner media border of image no. 16 reveals plaque, but is not delineated
by our algorithm. This lamellar plaque area can only be seen on the T2 image but cannot be
seen on the histological image. Subsequently, our algorithm does not find plaque.
5.2 Results Regarding the Healthy Media Border in the Local Weighting Concept
The Results for introducing the concept of the healthy media border in our algorithm
(according to equation (15)) are shown in Fig. 10.
18
Fig. 10: Results of the level set propagation of the local weighting concept supported by the knowledge of the
maximal shrinking distance according to equation (15). (The maximal shrinking distance is shown in magenta in
images of column 1 and 3).
19
The maximal shrinking distance avoids the leakage of the zero level set into plaque regions,
which are not identified by the deviant intensity areas and furthermore not delineated by
Canny edges (Please compare the results for image numbers 11-13 and 19 in Fig. 8, 9 and 10).
Obviously, the healthy media border serves as helpful information to improve the
segmentation results especially for vessel wall abnormalities which is not demarcated by
different signal levels in the MR images and thus cannot be recognized by a human observer
(Please note the result for image no. 15 in Fig. 10).
Such local alterations of vessel wall occur frequently in humans and are not yet regarded as
plaques, but described as positive remodeling (Glagov-effect). Those regions represent areas,
where plaque accumulation begins as expected and plaque tissue is most likely. For MR
images without plaque occurrence (i.e. image no. 18) the application of the maximal shrinking
distance could show plaque, which indeed represents vessel wall abnormalities. Subsequently,
in images without plaque the healthy media border does not detect plaque tissue, but reveals
vessel wall abnormalities, where plaque is most likely.
5.3 Comparison of the Results
Due to the difficult acquisition of this very exceptional MR-images we could only test our
algorithm for 12 different specimens, which we show in Fig. 5. Fortunately, the specimens
were measured in 256 slices. To avoid stronger correlations between those image data, we
selected 3 slices from every specimen, which are as far away as possible from each other.
Hence, we could test our level set method for 36 images and compare their results in this
chapter.
For analyzing the segmentation result we compared the segmented regions between the inner
media border and the final zero level set of the third level set step with the expert knowledge
about plaque for every image. We relate the sum of the incorrectly segmented areas to the
arterial area which represent the relative segmentation errors as shown in Fig. 11. For every
specimen number we clarify the segmentation errors with bars for the 3 selected slices. The
images which can be seen in Fig.5 are the most representative, hence their results are shown
with the 2 bars on the left for every specimen number. The results for the other slices can be
seen in the following 4 bars for every specimen number.
The results show that the most segmentation errors for the local weighting concept are lower
than 10 % except for the specimen numbers 15 and 19. However, if we compare the results
for the both applied methods, we see that the knowledge of the MSD does not reduce the
segmentation errors in general. Particularly for the specimen numbers 9, 10, 16, 17 and 18,
that do not show plaque or show only marginal plaque in all slices, the MSD algorithm
augments the segmentation errors definitely. For all the images coming from those specimen
numbers, the MSD algorithm reveals vessel wall abnormalities, which are recognized as
plaque tissue and increases misleadingly the segmentation errors.
The MSD algorithm reduces the segmentation errors for images that show higher plaque
complexity. Images that show more variety in neointimal plaque formations and chronic
plaque (specimen numbers 12, 14, 15 and 19) miss local deviant areas and edges. In those
20
cases the MSD can prevent the zero level set leaking into plaque regions. On the third slice of
specimen number 15, the vessel has separated and shows a completely different profile
For the application of the MSD, an exact demarcation of the healthy media border is required.
It is not determinable, if during preparation the vessel wall thickness has been altered locally.
This influences the delineation quality of the healthy media border. Generally, we recommend
the application of the local weighting concept with and without the MSD. Especially for
images without or with marginal plaque, the application of the algorithm without the MSD
controls, if the algorithm with MSD just detected vessel wall abnormalities. Furthermore, this
strategy delivers the extent of vessel wall abnormalities for images that do not show plaque
tissue.
For evaluating the segmentation results it should be mentioned, that the denoising technique
influences the segmentation results by removing thin plaque tissue structures which
consequently cannot be delineated by the level set algorithm applied to the denoised images.
Furthermore, the segmentation results are influenced by the indistinct media areas which seem
to be plaque and have been recognized by our algorithm, but cannot be determined
unambiguously as plaque in the histological images.
21
Fig. 11: Segmentation errors with respect to the area of the arterial profiles for the analyzed images. The two
applied methods are compared (lightgray: local weighting concept; black: local weighting concept with MSD).
22
6. Conclusion
We explain the preparation and measurement of the MR images we are working with in more
detail in Chapter 2. In Chapter 3 we describe the denoising technique we use to denoise the
MR images and prepare them for the level set propagation. At the beginning of Chapter 4 we
give theoretical backgrounds of the level set method. In the main part of this chapter we
present our method for the threestepped level set segmentation. For the first level set step we
use the geodesic active contour model and determine an initial seed point inside the lumen
(which should be determined automatically by applying the Hough transformation for circle
detection in future work). The geodesic active contour propagation guides the first extending
zero level set towards the exact boundaries of the inner media border. Due to strong intensity
gradients from plaque inside the media we adapt a level set technique [34] for Canny edges
instead of using the geodesic active contour concept. Basing on the final demarcation from the
first level set step we let the zero level set from the second step converge towards the outer
media border circulating around the thin Canny edges from plaque tissue. To segment the
media, we apply the level set technique on averaged and denoised T1 and T1fatsat images,
which exhibit sufficient edge magnitudes for approaching the inner and outer media border
with the zero level set.
The third level set step is based on the knowledge of the already segmented inner and outer
media border which we collected during the first and second level set segmentation phases.
By shrinking the zero level set, which is initially identical with the determined outer media
border, we approach towards the Canny edges of the denoised average image derived of
plaque relevant PD and T2 channels. Additionally, we propose a method to emphasize areas
with higher locally weighted intensity deviances to detect conspicuous plaque patterns. This
local weighting approach serves as a filter for plaque intensities and is used as a further
stopping criterion for the third zero level set. The local weighting concept evades the problem
of intensity attenuation along arbitrary directions in the image, which arises during the MR
measurement of the analyzed MR images. Furthermore, we are able to emphasize local
distinctive features by adjusting a new parameter for the standard deviation of the local
weighting function.
Finally, we discuss and justify the introduction of the maximal shrinking distance (MSD) and
incorporate this knowledge into the level set equation. With the MSD we can use a further
stopping criterion for the shrinking zero level set. The MSD delineates a region where plaque
is most likely and can improve the segmentation results.
We apply our algorithm for 36 images derived from 12 different specimens. The results of the
local weighting concept are shown in Chapter 5, where we show the segmentation results for
the threestepped algorithm with and without the MSD concept for the analyzed images. We
introduced expert knowledge and discussed which kind of expert knowledge is appropriate to
evaluate the segmentation results. In the same Chapter we list the segmentation errors based
upon the expert knowledge for all analyzed 36 images.
For the most images the segmentation errors for the local weighting concept are lower than 10
%. The results reveal, however, that for images with marginal or missing plaque the MSD
concept increases the segmentation errors. Here the vessel wall abnormalities are responsible
for this error augmentation. For images that show higher plaque complexities, the MSD
concept reduces the segmentation errors and serves as a reliable stopping criterion for the
third zero level set. We discuss representative results and give reasons for them. Furthermore,
23
we reason our recommendations for using both applied algorithms for the third step and give
hints how to interpret the results.
The application of the threestepped level set segmentation method shows good results,
although the available MR images are highly noisy and have been influenced by certain
preparations before their measurement in a modern MR tomograph. Due to these
circumstances we think, that adapted level set strategies for plaque recognition will also yield
useful segmentation results for worse in-plane resolutions from upcoming in-vivo plaque
measurements. We aim to encourage for the application of level set methods that use special
knowledge during the propagation (e.g. the MSD concept).
Upcoming, we aim to test iteratively the segmentation techniques for a successive
deterioration of in-plane resolutions to determine frontiers for image qualities in automatic
image understanding of plaque. Those contributions should improve the adjustment between
the research domains of automatic image understanding and biomedical MR-imaging. We
concentrate on level set segmentation with prior shape in future, when more images of similar
types of ex-vivo plaque images (which we analyzed in this work) are available. Furthermore,
we aim to improve the denoising methods with vector-valued denoising techniques to
preserve more of the thin plaque structures, which can be removed unintentionally by using
the scalar-valued anisotropic diffusion technique. To preserve all thin plaque structures we
applied the level set technique on the original noisy images. Unfortunately, the level set
propagation stopped too early due to high image gradients caused by low SNR. Although
those thin plaque patterns are actually of minor relevance we aim to enhance the denoising
results, which serves as an important preprocessing step for the application of the level set
segmentation and influences the segmentation results.
6. Acknowledgements
Special acknowledgements are directed to the members of the Department of Medical
Computer Science of Prof. Tolxdorff (FU Berlin) for all the helpful discussions. Especially,
we would like to thank Dr. Jürgen Braun for his contributions during the discussions. I would
like to thank Prof. Vicent Caselles for his very helpful explanations in the discussion we had
about level sets during SIMUMAT summer school 2007.
24
Literature
[1] Liu F, Xu D, Ferguson MS, Chu B, Saam T, Takaya N, Hatsukami TS, Yuan C, Kerwin
WS. Automated in vivo segmentation of carotid plaque MRI with Morphology-Enhanced
probability maps. Magnetic Resonance in Medicine 2005, 55(3), pp. 659- 668.
[2] Kim WY, Astrup AS, Stuber S, Tarnow L, Falk E, Botnar RM, Simonsen C, Pietraszek L,
Hansen PR, Manning WJ, Andersen NT, Parving HH. Subclinical coronary and aortic
atherosclerosis detected by magnetic resonance imaging in type 1 diabetes with and without
diabetic nephropathy. Journal of the American Heart Association 2007; 115(2): 228-35,
[3] Spuentrup E, Katoh M, Stuber M, Botnar RM, Schäffter T, Bücker A, Günther RW.
Coronary MR imaging using free-breathing 3D steady-state free precession with radial k-
space sampling. Fortschr Röntgenstr. 2003; 175, pp. 1330-4
[4] Spuentrup E , Ruebben A, Mahnken A, Stuber M, Kölker C, Trung HN, Günther RW,
Buecker A. Artifact-free coronary magnetic resonance angiography and coronary vessel wall
imaging in the presence of a new, metallic, coronary magnetic resonance imaging stent.
Journal of the American Heart Association 2005; 111, pp. 1019-1026.
[5] Jahnke C, Nehrke K, Paetsch I, Schnackenburg B, Gebker R, Fleck E, Nagel E. Improved
bulk myocardial motion suppression for navigator-gated coronary magnetic resonance
imaging, Journal Magnetic Resonance Imaging 2007, 26 (3), pp. 780-6.
[6] Langreck H, Schnackenburg B, Nehrke K, Boernert P, Wahl A, Paetsch I, Bornstedt A,
Fleck E, Nagel E. MR coronary artery imaging with 3D motion adapted gating (MAG) in
comparison to a standard prospective navigator technique. Journal Cardiovascular Magnetic
Resonance 2005; 7(5), pp. 793-7.
[7] Briley-Saebo KC, Mulder WJ, Mani V, Hyafil F, Amirbekian V, Aguinaldo JG, Fisher
EA, Fayad ZA. Magnetic resonance imaging of vulnerable atherosclerotic plaques: current
imaging strategies and molecular imaging probes. Journal Magnetic Resonance Imaging
2007; 26(3), 460-79.
[8] Frias JC, Williams KJ, Fisher EA, Fayad ZA. Recombinant HDL-like nanoparticles: a
specific contrast agent for MRI of atherosclerotic plaques. Journal American Chemical
Society 2004; 126(50), 16316-16317.
[9] Gloger O, Jaeger M, Nagel E, Hellwich O. Hierarchical Scale Space Pre-Segmentation
based on Attraction Areas for MR-images from Atherosclerosis-plaque. Proc. of International
Workshop on Image Analysis in the Life Sciences Theory and Applications (Linz) 2007, p.22.
[10] Kass M, Witkin A, Terzopoulos D. Snakes: Active contour models. International
Journal of Computer Vision 1988, 1, pp.321-331.
[11] Cohen LD, Cohen I. Finite element methods for active contour models and balloons for
2D and 3D images. IEEE Trans. on Pattern Analysis and Machine Intelligence 1993; 15(11),
1131–1147.
[12] Davatzikos CA, Prince JL, An Active Contour Model for Mapping the Cortex. IEEE
Transactions on Medical Imaging 1995; 14(1), 65-80.
25
[13] Desbleds Mansard C, Canet Soulas EP, Anwander A, Chaabane L, Neyran B, Serfaty
JM, Magnin IE, Douek PC, Orkisz M. Quantification of multi-contrast vascular MR Images
with the NLSnake, an active contour model: in vitro validation and in vivo evaluation.
Magnetic Resonance in Medicine 2004; 51(2), 370-379.
[14] Rougon N, Prˆeteux F. Directional adaptive deformable models for segmentation with
application to 2D and 3D medical images. Medical Imaging 93: Image Processing 1993; 1898
of SPIE Proc., 193–207. Bellingham, WA: SPIE.
[15] Singh A, von Kurowski L, Chiu MY. Cardiac MR image segmentation using deformable
models. Biomedical Image Processing and Biomedical Visualization 1993; 1905 of SPIE
Proc., 8–28. Bellingham, WA: SPIE.
[16] Xu C, Prince JL. Gradient Vector Flow: A New External Force for Snakes. IEEE
Proceeding Conference On Computer Vision and Pattern Recognition 1997; pp. 66-71.
[17] Caselles V, Kimmel R, Sapiro G. Geodesic Active Contours. International Journal of
Computer Vision 1997; 22(1), 61-79.
[18] Rousson M, Deriche R. A Variational Framework for Active and Adaptive Segmentation
of Vector Valued Images. In Proc. IEEE Workshop on Motion and Video Computing,
Orlando, Florida, Dec. 2002.
[19] Chan T, Vese L. “Active Contours without Edges”, IEEE Transactions on Image
Processing, Vol. 10 (2), 2001, pp. 266-277.
[20] Chan T, Vese L. A MultiPhase Level Set Framework for Image Segmentation Using the
Mumford and Shah Model. International Journal of Computer Vision, Vol. V50, No. 3. (1
December 2002), pp. 271-293.
[21] Sum KW, Paul YS. A Novel Active Contour Model Using Local and Global Statistics
for Vessel Extraction. IEEE Engineering in Medicine and Biology Society, 2006. EMBS 06,
28th Annual International Conference of the IEEE, pp. 3126-3129.
[22] Lankton S, Nain D, Yezzi A, Tannenbaum A. Hybrid geodesic region-based curve
evolutions for image segmentation. Medical Imaging 2007: Physics of Medical Imaging.
Edited by Hsieh, Jiang; Flynn, Michael J. Proceedings of the SPIE, Volume 6510, pp. 65104U
(2007).
[23] Jahnke C, Dietrich T, Paetsch I, Koehler U, Preetz K, Schnackenburg B, Fleck E, Graf K,
Nagel E. Experimental evaluation of the detectability of submillimeter atherosclerotic lesions
in ex vivo human iliac arteries with ultrahigh-field (7.0 T) magnetic resonance imaging.
International Journal for Cardiovascular Imaging 2006; 23(4), 519-527
[24] Perona P, Malik J, Scale space and edge detection using anisotropic diffusion, IEEE
Trans. Pattern Anal. Mach. Intell. 1990; 12, pp. 629-639.
[25] Wei GW. Generalized Perona-Malik equation for Image Restoration, IEEE Signal
Processing Lett. , Vol 6, 165-167 (1999)
26
27
[26] Whitaker RT, Xue X. Variable-Conductance, Level-Set Curvature for Image Processing,
IEEE International Conference on Image Processing 2001; pp. 142–145.
[27] Sapiro G. Vector (self) snakes: a geometric framework for color, texture, and multiscale
image segmentation. Proceedings IEEE Int. Conf. Image Processing (ICIP-96, Lausanne,
Sept. 16-19, 1996), 1, pp. 817-820.
[28] Ibanez L, Schroeder W. The ITK Software Guide, Kitware, Inc., 2005.
[29] Dervieux A, Thomasset F. A finite element method for the simulation of a Rayleigh-
Taylor instability, In: R.Rautman (Ed.): Approximation Methods for Navier-Stokes Problems,
Lecture Notes in Mathematics, 771, 145-158. Springer, Berlin, 1979.
[30] Osher S, Sethain JA. Fronts propagating with curvature dependent speed: Algorithms
based on Hamilton-Jacobi formulation. Journal of Computational Physics 1988; 79, 12-49.
[31] Osher S, Fedkiw R. Level Set Methods and Dynamic Implicit Surfaces, Springer-Verlag,
2003
[32] Osher S, Paragios N. Geometric Level Set Methods in Imaging, Vision and Graphics,
Springer-Verlag New York, 2003
[33] Sethian JA. Level Set Methods and Fast Marching Methods, Cambridge University
Press, 1999
[34] Malladi R, Sethian J, Vemuri B. Shape modelling with front propagation: A level set
approach. IEEE Trans. on Pattern Analysis and Machine Intelligence 1995; 17(2), 158-174.
[35] Caselles V, Catte F, Coll T, Dibos F. A geometric model for active contours in image
processing. Numerische Mathematik 1993; 66, 1-31.
[36] Canny J, A computational approach to edge detection. IEEE Trans. on Pattern Analysis
and Machine Intelligence 1986; 8(6), 679-698.
[37] Qin X, Jiang J, Wang W, Zhang F. Canny Operator Based Level Set Segmentation
Algorithm for Medical Images, Bioinformatics and Biomedical Engineering 2007. First
International Conference on Volume, 6-8 July 2007 pp. 892 – 895.
