ABSTRACT
Detection of Regions of Interest
Although a physician or a radiologist of necessity will carefully examine
an image on hand in its entirety more often than not diagnostic features of
interest manifest themselves in local regions It is uncommon that a condition
or disease will alter an image over its entire spatial extent In a screening
situation the radiologist scans the entire image and searches for features
that could be associated with disease In a diagnostic situation the medical
expert concentrates on the region of suspected abnormality and examines its
characteristics to decide if the region exhibits signs related to a particular
disease
In the CAD environment one of the roles of image processing would be to
detect the region of interest ROI for a given specic screening or diagnostic
application Once the ROIs have been detected the subsequent tasks would
relate to the characterization of the regions and their classication into one
of several categories A few examples of ROIs in dierent biomedical imaging
and image analysis applications are listed below
Cells in cervicalsmear test images Papanicolaou or Papsmear test
Calcications in mammograms
Tumors and masses in mammograms
The pectoral muscle in mammograms
The breast outline or skinair boundary in mammograms
The broglandular disc in mammograms
The airway tree in lungs
The arterial tree in lungs
The arterial tree of the left ventricle and constricted parts of the same
due to plaque development
Segmentation is the process that divides an image into its constituent parts
objects or ROIs Segmentation is an essential step before the description
recognition or classication of an image or its constituents Two major ap
proaches to image segmentation are based on the detection of the following
characteristics
DiscontinuityAbrupt changes in gray level corresponding to edges
are detected
Similarity Homogeneous parts are detected based on graylevel
thresholding region growing and region splittingmerging
Depending upon the nature of the images and the ROIs we may attempt to
detect the edges of the ROIs if distinct edges are present or we may attempt
to grow regions to approximate the ROIs It should be borne in mind that
in some cases an ROI may be composed of several disjoint component areas
for example a tumor that has metastasized into neighboring regions and
calcications in a cluster Edges that are detected may include disconnected
parts that may have to be matched and joined We shall explore several
techniques of this nature in the present chapter
Notwithstanding the stated interest in local regions as above applications
do exist where entire images need to be analyzed for global changes in pat
terns for example changes in the orientational structure of collagen bers in
ligaments see Figure and bilateral asymmetry in mammograms see Sec
tion Furthermore in the case of clustered calcications in mammograms
cells in cervical smears and other examples of images with multicomponent
ROIs analysis may commence with the detection of single units of the pattern
of interest but several such units present in a given image may need to be
analyzed separately and together in order to reach a decision regarding the
case
Thresholding and Binarization
If the gray levels of the objects of interest in an image are known from prior
knowledge or can be determined from the histogram of the given image the
image may be thresholded to detect the features of interest and reject other
details For example if it is known that the objects of interest in the image
have graylevel values greater than L
we could create a binary image for
display as
gmn
if fmn L
if fmn L
where fmn is the original image gmn is the thresholded image to be
displayed and the display range is See also Section
Methods for the derivation of optimal thresholds are described in Sections
and
Example Figure a shows a TEM image of a ligament sample demon
strating collagen bers in crosssection see Section Inspection of the his
togram of the image shown in Figure shows that the sections of the
collagen bers in the image have graylevel values less than about values
greater than this level represent the brighter background in the image The
histogram also indicates the fact that the graylevel ranges of the collagen
ber regions and the background overlap signicantly Figure b shows a
thresholded version of the image in a with all pixels less than appearing
in black and all pixels above this level appearing in white This operation
is the same as the thresholding operation given by Equation but in the
opposite sense Most of the collagen ber sections have been detected by the
thresholding operation However some of the segmented regions are incom
plete or contain holes whereas some parts that appear to be separate and
distinct in the original image have been merged in the result An optimal
threshold derived using the methods described in Sections and
could lead to better results
Detection of Isolated Points and Lines
Isolated points may exist in images due to noise or due to the presence of
small particles in the image The detection of isolated points is useful in noise
removal and the analysis of particles The following convolution mask may
be used to detect isolated points
The operation computes the dierence between the current pixel at the center
of the mask and the average of its connected neighbors The mask could
also be seen as a generalized version of the Laplacian mask in Equation
The result of the mask operation could be thresholded to detect isolated pixels
where the dierence computed would be large
Straight lines or line segments oriented at
o
o
o
and
o
may be
detected by using the following convolution masks
A line may be said to exist in the direction for which the corresponding mask
provides the largest response
a
b
FIGURE
a TEM image of collagen bers in a scartissue sample from a rabbit ligament
at a magnication of approximately See also Figure Image
courtesy of CB Frank Department of Surgery University of Calgary See
Figure for the histogram of the image b Image in a thresholded at
the gray level of
Edge Detection
One of the approaches to the detection of an ROI is to detect its edges The
HVS is particularly sensitive to edges and gradients and some theories and
experiments indicate that the detection of edges plays an important role in
the detection of objects and analysis of scenes
In Section on the properties of the Fourier transform we saw that
the rstorder derivatives and the Laplacian relate to the edges in the im
age Furthermore we saw that these spacedomain operators have equivalent
formulations in the frequency domain as highpass lters with gain that is
proportional to frequency in a linear or quadratic manner The enhancement
techniques described in Sections and further strengthen the relation
ship between edges gradients and highfrequency spectral components We
shall now explore how these approaches may be extended to detect the edges
or contours of objects or regions
Note Some authors consider edge extraction to be a type of image en
hancement
Convolution mask operators for edge detection
An edge is characterized by a large change in the gray level from one side
to the other in a particular direction dependent upon the orientation of the
edge Gradients or derivatives measure the rate of change and hence could
serve as the basis for the development of methods for edge detection
The rst derivatives in the x and y directions approximated by the rst
dierences are given by using matrix notation
f
yb
mn fmn fm n
f
xb
mn fmn fmn
where the additional subscript b indicates a backwarddierence operation
Because causality is usually not a matter of concern in image processing the
dierences may also be dened as
f
yf
mn fm n fmn
f
xf
mn fmn fmn
where the additional subscript f indicates a forwarddierence operation A
limitation of the operators as above is that they are based upon the values
of only two pixels this makes the operators susceptible to noise or spurious
pixel values A simple approach to design robust operators and reduce the
sensitivity to noise is to incorporate averaging over multiple measurements
Averaging the two denitions of the derivatives in Equations and we
get
f
ya
mn fm n fm n
f
xa
mn fmn fmn
where the additional subscript a indicates the inclusion of averaging
In image processing it is also desirable to express operators in terms of
oddsized masks that may be centered upon the pixel being processed The
Prewitt operators take these considerations into account with the following
masks for the horizontal and vertical derivatives G
x
and G
y
respectively
G
x
G
y
The Prewitt operators use three dierences across pairs of pixels in three
rows or columns around the pixel being processed Due to this fact and due
to the scale factor of in Equation in order to derive the exact gradient
the results of the Prewitt operators should be divided by where
is the sampling interval in x and y however this step could be ignored if the
result is scaled for display or thresholded to detect edges
In order to accommodate the orientation of the edge a vectorial form of
the gradient could be composed as
G
f
mn G
fx
mn j G
fy
mn
kG
f
mnk
q
G
fx
mn G
fy
mn
G
f
mn tan
G
fy
mn
G
fx
mn
where
G
fx
mn f G
x
mn
and
G
fy
mn f G
y
mn
If the magnitude is to be scaled for display or thresholded for the detec
tion of edges the squareroot operation may be dropped or the magnitude
approximated as jG
fx
j jG
fy
j in order to save computation
The Sobel operators are similar to the Prewitt operators but include larger
weights for the pixels in the row or column of the pixel being processed as
G
x
G
y
Edges oriented at
o
and
o
may be detected by using rotated versions
of the masks as above The Prewitt operators for the detection of diagonal
edges are
G
o
and
G
o
Similar masks may be derived for the Sobel operator
Note The positive and negative signs of the elements in the masks above
may be interchanged to obtain operators that detect gradients in the opposite
directions This step is not necessary if directions are considered in the range
o
o
only or if only the magnitudes of the gradients are required
Observe that the sum of all of the weights in the masks above is zero
This indicates that the operation being performed is a derivative or gradient
operation which leads to zero output values in areas of constant gray level
and the loss of intensity information
The Roberts operator uses neighborhoods to compute crossdierences
as
and
The masks are positioned with the upperleft element placed on the pixel
being processed The absolute values of the results of the two operators are
added to obtain the net gradient
gmn jfm n fmnj jfm n fmn j
with the indices in matrixindexing notation The individual dierences may
also be squared and the square root of their sum taken to be the net gradient
The advantage of the Roberts operator is that it is a forwardlooking operator
as a result of which the result may be written in the same array as the input
image This was advantageous when computer memory was expensive and in
short supply
Examples Figure a shows the Shapes text image Part b of the
gure shows the gradient magnitude image obtained by combining as in
Equation the horizontal and vertical derivatives shown in parts c and
d of the gure respectively The image in part c presents high values pos
itive or negative at vertical edges only horizontally oriented edges have been
deleted by the horizontal derivative operator The image in part d shows
high output at horizontal edges with the vertically oriented edges having been
removed by the vertical derivative operator The test image has strong edges
for most of the objects present which are clearly depicted in the derivative
images however the derivative images show the edges of a few objects that
are not readily apparent in the original image as well Parts e and f of the
gure show the derivatives at
o
and
o
respectively the images indicate
the diagonal edges present in the image
Figures and show similar sets of results for the clock the
knee MR and the chest Xray test images respectively In the derivatives
of the clock image observe that the numeral has been obliterated by
the vertical derivative operator Figure d but gives rise to high output
values for the horizontal derivative Figure c The clock image has the
minute hand oriented at approximately
o
with respect to the horizontal
this feature has been completely removed by the
o
derivative operator as
shown in Figure f but has been enhanced by the
o
derivative operator
as shown in Figure e The knee MR image contains sharp boundaries
that are depicted well in the derivative images in Figure The derivative
images of the chest Xray image in Figure indicate large values at the
boundaries of the image but depict the internal details with weak derivative
values indicative of the smooth nature of the image
The Laplacian of Gaussian
Although the Laplacian is a gradient operator it should be recognized that it
is a secondorder dierence operator As we observed in Sections and
this leads to doubleedged outputs with positive and negative values at
each edge this property is demonstrated further by the example in Figure
see also Figure The Laplacian has the advantage of being omnidirec
tional that is being sensitive to edges in all directions however it is not
possible to derive the angle of an edge from the result The operator is also
sensitive to noise because there is no averaging included in the operator the
gain in the frequency domain increases quadratically with frequency caus
ing signicant amplication of highfrequency noise components For these
reasons the Laplacian is not directly useful in edge detection
The doubleedged output of the Laplacian indicates an important property
of the operator the result possesses a zerocrossing in between the positive
and negative outputs across an edge the property holds even when the edge
in the original image is signicantly blurred This property is useful in the
development of robust edge detectors The noise sensitivity of the Laplacian
a b
c d
e f
FIGURE
a Shapes test image b Gradient magnitude display range out of
c Horizontal derivative display range out of
d Vertical derivative display range out of e
o
derivative display range out of f
o
derivative
display range out of
a b
c d
e f
FIGURE
a Clock test image b Gradient magnitude display range out of
c Horizontal derivative display range out of
d Vertical derivative display range out of e
o
derivative display range out of f
o
derivative
display range out of
a b
c d
e f
FIGURE
a Knee MR image b Gradient magnitude display range out of
c Horizontal derivative display range out of
d Vertical derivative display range out of e
o
derivative display range out of f
o
derivative
display range out of
a b
c d
e f
FIGURE
a Part of a chest Xray image b Gradient magnitude display range
out of c Horizontal derivative display range out of
d Vertical derivative display range out of
e
o
derivative display range out of f
o
deriva
tive display range out of
FIGURE
Top to bottom A prole of a blurred object showing two edges the rst
derivative and the second derivative see also Figure
may be reduced by including a smoothing operator A scalable smoothing op
erator could be dened in terms of a D Gaussian function with the variance
controlling the spatial extent or width of the smoothing function Combining
the Laplacian and the Gaussian we obtain the popular LaplacianofGaussian
or LoG operator
Consider the Gaussian specied by the function
gx y exp
x
y
The usual normalizing scale factor has been left out Taking partial derivatives
with respect to x and y we obtain
g
x
x
exp
x
y
g
y
y
exp
x
y
which leads to
r
gx y LoGr
r
exp
r
where r
p
x
y
The LoG function is isotropic and has positive and
negative values Due to its shape it is often referred to as the Mexican hat
or sombrero
Figure shows the LoG operator in image and meshplot formats the
basic Gaussian used to derive the LoG function is also shown for reference
The Fourier magnitude spectra of the Gaussian and LoG functions are also
shown in the gure It should be observed that whereas the Gaussian is a
lowpass lter which is also a D Gaussian in the frequency domain the LoG
function is a bandpass lter The width of the lters is controlled by the
parameter of the Gaussian
Figure shows proles of the LoG and the related Gaussian for two values
of Figure shows the proles of the Fourier transforms of the functions
in Figure The proles clearly demonstrate the nature of the functions
and their ltering characteristics
An approximation to the LoG operator is provided by taking the dierence
between two Gaussians of appropriate variances this operator is known as
the dierenceofGaussians or DoG operator
Examples Figure shows the Shapes test image the LoG of the image
with pixel and the locations of the zerocrossings of the LoG of the
image with pixel and pixels The zerocrossings indicate the
locations of the edges in the image The use of a large value for reduces the
eect of noise but also causes smoothing of the edges and corners as well as
the loss of the minor details present in the image
a b
c d
e f
FIGURE
The Laplacian of Gaussian in b image format and d as a mesh plot The
related Gaussian functions are shown in a and c The size of the arrays
is pixels standard deviation pixels The Fourier magnitude
spectra of the functions are shown in e and f
a
b
FIGURE
Proles of the Laplacian of Gaussian solid line and the related Gaussian
dashed line in Figure The functions have been normalized to a maxi
mum value of unity The unit of r is pixels a pixels b pixels
a
b
FIGURE
Proles of the Fourier magnitude spectra of the Laplacian of Gaussian solid
line and the related Gaussian dashed line in Figure Both functions
have been normalized to have a maximum value equal to unity a
pixels b pixels The zerofrequency point is at the center of the
horizontal axis
Figure shows the clock image its LoG and the zerocrossings of the
LoG with pixel and pixels The results illustrate the performance
of the LoG operator in the presence of noise
Figures and show similar sets of results for the myocyte
image the knee MR image and the chest Xray test images Comparative
analysis of the scales of the details present in the images and the zerocrossings
of the LoG for dierent values of indicates the importance of selecting values
of the parameter in accordance with the scale of the details to be detected
Scalespace methods for multiscale edge detection
Marr and Hildreth suggested that physical phenomena may be de
tected simultaneously over several channels tuned to dierent spatial sizes
or scales with an approach known as the spatial coincidence An intensity
change that is due to a single physical phenomenon is indicated by zero
crossing segments present in independent channels over a certain range of
scales with the segments having the same position and orientation in each
channel A signicant intensity change indicates the presence of a major event
that is registered as a physical boundary and is recognized as a single phys
ical phenomenon The boundaries of a signicant physical pattern should be
present over several channels suggesting that the use of techniques based on
zerocrossings generated from lters of dierent scales could be more eec
tive than the conventional singlescale methods for edge detection see for
example Figure
Zerocrossings and scalespace The multichannel model for the HVS
and the MarrHildreth spatial coincidence assumption led to the
development of methods for the detection of edges based upon multiscale
analysis performed with lters of dierent scales Marr and Hildreth pro
posed heuristic rules to combine information from the dierent channels in a
multichannel vision model they suggested the use of a bank of LoG lters
with several values of which may be represented as fr
gx yg with
A method for obtaining information in images across a continuum of scales
was suggested byWitkin who introduced the concept of scalespace The
method rapidly gained considerable interest and has been explored further by
several researchers in image processing and analysis
The scalespace x y of an image fx y is dened as the set of all zero
crossings of its LoG
fx yg fx yg j x y
x
y
where
x y fr
gx y fx yg
a b
c d
FIGURE
a The Shapes test image b The LoG of the image in a with pixel
c Locations of the zerocrossings in the LoG in b d Locations of the
zerocrossings in the LoG with pixels
a b
c d
FIGURE
a The clock test image b The LoG of the image in a with pixel
c Locations of the zerocrossings in the LoG in b d Locations of the
zerocrossings in the LoG with pixels
a b
c d
FIGURE
a Image of a myocyte b The LoG of the image in a with pixels
c Locations of the zerocrossings in the LoG in b d Locations of the
zerocrossings in the LoG with pixels
a b
c d
e f
FIGURE
a MR image of a knee b The LoG of the image in a with pixels
c f Locations of the zerocrossings in the LoG with and
pixels
a b
c d
FIGURE
a Part of a chest Xray image b The LoG of the image in a with
pixels display range out of c Locations of the zero
crossings in the LoG in b d Locations of the zerocrossings in the LoG
with pixels
As the scale varies from to the set fx yg forms continuous
surfaces in the x y scalespace
Several important scalespace concepts apply to D and D signals It has
been shown that the scalespace of almost all signals ltered by a Gaussian
determines the signal uniquely up to a scaling constant except for noise
contaminated signals and some special functions The importance of
this property lies in the fact that theoretically for almost all signals no
information is lost by working in the scalespace instead of the image domain
This property plays an important role in image understanding image
reconstruction from zerocrossings and image analysis using the
scalespace approach Furthermore it has also been shown that the
Gaussian does not create additional zerocrossings as the scale increases
beyond a certain limit and that the Gaussian is the only lter with this
desirable scaling behavior
Based on the spatialcoincidence assumption Witkin proposed a D
stability analysis method for the extraction of primitive events that occur over
a large range of scales The primitive events were organized into a qualitative
signal description representing the major events in the signal Assuming that
zerocrossing curves do not cross one another which was later proven to be
incorrect by Katz Witkin dened the stability of a signal interval as
the scale range over which the signal interval exists major events could then
be captured via stability analysis However due to the complex topological
nature of spatial zerocrossings it is often dicult to directly extend Witkins
D stability analysis method to D image analysis The following problems
aect Witkins method for stability analysis
It has been shown that zerocrossing curves do cross one another
It has been shown that real authentic zerocrossings could turn into
false phantom zerocrossings as the scale increases Use of
the complete scalespace with ranging from to may introduce
errors in certain applications an appropriate scalespace using only a
nite range of scales could be more eective
For D signals images the scalespace consists of zerocrossing surfaces
that are more complex than the zerocrossing curves for D signals The
zerocrossing surfaces may split and merge as the scale varies decreases
or increases respectively
As a consequence of the above there is no simple topological region as
sociated with a zerocrossing surface and tracing a zerocrossing surface
across scales becomes computationally dicult
Liu et al proposed an alternative denition of zerocrossing surface
stability in terms of important spatial boundaries In this approach a spatial
boundary is dened as a region of steep gradient and high contrast and is
welldened if it has no neighboring boundaries within a given range This
denition of spatial boundaries is consistent with the MarrHildreth spatial
coincidence assumption Furthermore stability maps associated with
the scalespace are used A relaxation algorithm is included in the process to
generate zerocrossing maps
In the method of Liu et al the discrete scalespace approach is used to
construct a representation of a given image in terms of a stability map which
is a measure of pattern boundary persistence over a range of lter scales For
a given image fx y a set of zerocrossing maps is generated by convolving
the image with the set of isotropic functions r
gx y
i
i N It
was indicated that N sampled
i
values ranging from to pixels were
adequate for the application considered Ideally one would expect a pattern
boundary to be accurately located over all of the scales However it has been
shown that the accuracy of zerocrossing localization depends upon
the width of the central excitatory region of the lter dened as w
i
p
i
Chen and Medioni proposed a D method for localization of
zerocrossings that works well for ideal step edges and image patterns with
sharp contrast however the method may not be eective for the construction
of the spatial scalespace for reallife images with poor and variable contrast
Instead of directly matching all the zerocrossing locations at a point x y
over the zerocrossing maps Liu et al proposed a criterion C
i
that is
a function of the scale
i
to dene a neighborhood in which the matching
procedure is performed at a particular scale
C
i
fx
y
g j x
i
x
x
i
y
i
y
y
i
where x
y
are the actual locations of the zerocrossings x y is the pixel
location at which the lters are being applied and is a constant to be
determined experimentally was used by Liu et al Therefore if a
zerocrossing x y
i
is found in the neighborhood dened by C
i
an
arbitrary constant is assigned to a function S
i
x y which otherwise is
assigned a zero that is
S
i
x y
if x y
i
C
i
otherwise
where the subscript i corresponds to the i
th
scale
i
Applying Equations and to the set of zerocrossings detected a set
of adjusted zerocrossing maps fS
x y S
x y S
N
x yg is obtained
where N is the number of scales The adjusted zerocrossing maps are used
to construct the zerocrossing stability map x y as
x y
N
X
i
S
i
x y
The values of x y are in principle a measure of boundary stability through
the lter scales Marr and Hildreth and Marr and Poggio suggested
that directional detection of zerocrossings be performed after the LoG oper
ator has been applied to the image
According to the spatialcoincidence assumption a true boundary should
be high in contrast and have relatively large values at the corresponding
locations Furthermore there should be no other edges within a given neigh
borhood Thus if in a neighborhood of x y nonzero stability map values
exist only along the orientation of a local segment of the stability map that
crosses x y then x y may be considered to signify a stable edge pixel at
x y On the other hand if many nonzero stability map values are present at
dierent directions x y indicates an insignicant boundary pixel at x y
In other words a consistent stability indexing method in the sense of the
spatialcoincidence assumption should take neighboring stability indices into
account Based upon this argument Liu et al proposed a relative stability
index x y computed from the stability map where x y as follows
In a neighborhood of x y if m nonzero values are found x y is
relabeled as l
and the rest of x
k
y
k
are relabeled as l
k
k m
see Figure In order to avoid using elements in the neighborhood
that belong to the same edge those x
y
having the same orientation as
that of l
are not included in the computation of x y Based upon these
requirements the relative stability index x y is dened as
x y
l
P
m
k
k
l
k
where
k
expd
k
and d
k
p
x x
k
y y
k
and x
k
y
k
are the
locations of l
k
It should be noted that x y and that the value of
x y is governed by the geometrical distribution of the neighboring stability
index values
Stability of zerocrossings Liu et al observed that the use of zero
crossings to indicate the presence of edges is reliable as long as the edges are
wellseparated otherwise the problem of false zerocrossings could arise
The problem with using zerocrossings to localize edges is that the zeros of
the second derivative of a function localize the extrema in the rst derivative
of the function the extrema include the local minima and maxima in the
rst derivative of the function whereas only the local maxima indicate the
presence of edge points Intuitively those zerocrossings that correspond to
the minima of the rst derivative are not associated with edge points at all
In image analysis based upon the notion of zerocrossings it is desirable to be
able to distinguish real zerocrossings from false ones and to discard the false
zerocrossings Motivated by the work of Richter and Ullman Clark
conducted an extensive study on the problem of false and real zerocrossings
and proposed that zerocrossings may be classied as real if x y and
false if x y where x y rr
px y rpx y where denotes
the dot product px y is a smoothed version of the given image such as
px y gx y fx y rpx y
h
p
x
p
y
i
T
and rr
px y
FIGURE
A case where three zerocrossings fl
l
l
g are found in a neighborhood of a
zerocrossing l
d
i
indicates the distance from l
i
to l
The arrows indicate
the directions of the zerocrossings a The neighboring zerocrossings are far
apart from l
imposing a low penalty to the zerocrossing associated with l
b The neighboring zerocrossings are close to l
imposing a high penalty to
the zerocrossing associated with l
Reproduced with permission from ZQ
Liu RM Rangayyan and CB Frank Statistical analysis of collagen align
ment in ligaments by scalespace analysis IEEE Transactions on Biomedical
Engineering
c
IEEE
h
p
x
p
x y
p
x
y
p
y
i
T
Liu et al included this step in their method
to detect true zerocrossings
Example Figure a shows an SEM image of collagen bers in a lig
ament scartissue sample Parts b d of the gure show the zerocrossing
maps obtained with and pixels The result in b contains several
spurious or insignicant zerocrossings whereas that in d contains smoothed
edges of only the major regions of the image Part e shows the stability map
which indicates the edges of the major objects in the image The stability map
was used to detect the collagen bers in the image See Section for details
on directional analysis of oriented patterns by further processing of the stabil
ity map Methods for the directional analysis of collagen bers are described
in Section
See Sections and for more examples on multiscale analysis
Cannys method for edge detection
Canny proposed an approach for edge detection based upon three criteria
for good edge detection multidirectional derivatives multiscale analysis and
optimization procedures The three criteria relate to low probabilities of false
edge detection and missing real edges represented in the form of an SNR good
localization represented by the RMS distance of the detected edge from the
true edge and the production of a single output for a single edge represented
by the distance between the adjacent maxima in the output A basic lter
derived using the criteria mentioned above was approximated by the rst
derivative of a Gaussian Procedures were proposed to incorporate multiscale
analysis and directional lters to facilitate ecient detection of edges at all
orientations and scales including adaptive thresholding with hysteresis
The LoG lter is nondirectional whereas Cannys method selectively evalu
ates a directional derivative across each edge By avoiding derivatives at other
angles that would not contribute to edge detection but increase the eects of
noise Cannys method could lead to better results than the LoG lter
Fourierdomain methods for edge detection
In Section we saw that highpass lters may be applied in the Fourier
domain to extract the edges in the given image However the inclusion of all
of the highfrequency components present in the image could lead to noisy
results Reduction of highfrequency noise suggests the use of bandpass l
ters which may be easily implemented as a cascade of a lowpass lter with a
highpass lter In the frequency domain such a cascade of lters results in
the multiplication of the corresponding transfer functions Because edges are
often weak or blurred in images some form of enhancement of the correspond
ing frequency components would also be desirable This argument leads us
to the LoG lter a combination of the Laplacian which is a highfrequency
FIGURE
a SEM image of a ligament scartissue sample b d Zerocrossing loca
tions detected using the LoG operator with and pixels respectively
e The stability map depicting the major edges present in the image Re
produced with permission from ZQ Liu RM Rangayyan and CB Frank
Statistical analysis of collagen alignment in ligaments by scalespace anal
ysis IEEE Transactions on Biomedical Engineering
c
IEEE
emphasis lter with its gain quadratically proportional to frequency with a
Gaussian lowpass lter The methods and results presented in Section
demonstrate the edgedetection capabilities of the LoG lter which may be
easily implemented in the frequency domain Frequencydomain implemen
tation using the FFT may be computationally advantageous when the LoG
function is specied with a large spatial array which would be required in the
case of large values of
Several other linedetection and edgedetection methods such as Gabor
lters see Sections and and fan lters see Section
may also be implemented in the frequency domain with advantages
Edge linking
The results of most methods for edge detection are almost always discontin
uous and need to be processed further to link disjoint segments and obtain
complete representations of the boundaries of ROIs Two principal proper
ties that may be used to establish the similarity of edge pixels from gradient
images are the following
The strength of the gradient a point x
y
in a neighborhood of
x y is similar in gradient magnitude to the point x y if
kGx yGx
y
k T
where Gx y is the gradient vector of the given image fx y at x y
and T is a threshold
The direction of the gradient a point x
y
in a neighborhood of
x y is similar in gradient direction to the point x y if
jx y x
y
j A
x y
Gx y tan
fx yy
fx yx
where A is a threshold
or neighborhoods may be used for checking pixels for similarity in
their gradients as above Further processing steps may include linking of edge
segments separated by small breaks and deleting isolated short segments
See Section page for the details of an edge analysis method
known as edge ow propagation
Segmentation and Region Growing
Dividing an image into regions that could correspond to structural units ob
jects of interest or ROIs is an important prerequisite for most techniques for
image analysis Whereas a human observer may by merely looking at a dis
played image readily recognize its structural components computer analysis
of an image requires algorithmic analysis of the array of image pixel values
before arriving at conclusions about the content of the image Computer
analysis of images usually starts with segmentation which reduces pixel data
to regionbased information about the objects and structures present in the
image
Image segmentation techniques may be classied into four main categories
thresholding techniques
boundarybased methods
regionbased methods and
hybrid techniques that combine boundary and
region criteria
Thresholding methods are based upon the assumption that all pixels whose
values lie within a certain range belong to the same class see Section
The threshold may be determined based upon the valleys in the histogram of
the image however identifying thresholds to segment objects is not easy even
with optimal thresholding techniques Moreover because threshold
ing algorithms are solely based upon pixel values and neglect all of the spatial
information in the image their accuracy of segmentation is limited further
more thresholding algorithms do not cope well with noise or blurring at object
boundaries
Boundarybased techniques make use of the property that usually pixel
values change rapidly at the boundaries between regions The methods start
by detecting intensity discontinuities lying at the boundaries between ob
jects and their backgrounds typically through a gradient operation High
values of the output provide candidate pixels for region boundaries which
must then be processed to produce closed curves representing the boundaries
between regions as well as to remove the eects of noise and discontinu
ities due to nonuniform illumination and other eects Although edgelinking
algorithms have been proposed to assemble edge pixels into a meaningful
set of object boundaries see Section such as local similarity analy
sis Houghtransformbased global analysis and global processing via graph
theoretic techniques the accurate conversion of disjoint sets of edge pixels
to closedloop boundaries of ROIs is a dicult task
Regionbased methods which are complements of the boundarybased ap
proach rely on the postulate that neighboring pixels within a region have
similar values Regionbased segmentation algorithms may be divided into
two groups region splitting and merging and region growing
Segmentation techniques in the region splitting and merging category ini
tially subdivide the given image into a set of arbitrary disjoint regions and
then merge andor split the regions in an attempt to satisfy some prespecied
conditions
Region growing is a procedure that groups pixels into regions The simplest
of regiongrowing approaches is pixel aggregation which starts with a seed
pixel and grows a region by appending spatially connected neighboring pixels
that meet a certain homogeneity criterion Dierent homogeneity criteria
will lead to regions with dierent characteristics It is important as well as
dicult to select an appropriate homogeneity criterion in order to obtain
regions that are appropriate for the application on hand
Typical algorithms in the group of hybrid techniques rene image segmen
tation by integration of boundary and region information proper combination
of boundary and region information may produce better segmentation results
than those obtained by either method on its own For example the mor
phological watershed method is generally applied to a gradient image
which can be viewed as a topographic map with boundaries between regions
as ridges Consequently segmentation is equivalent to ooding the topogra
phy from the seed pixels with region boundaries being erected to keep water
from the dierent seed pixels from merging Such an algorithm is guaranteed
to produce closed boundaries which is known to be a major problem with
boundarybased methods However because the success of this type of an al
gorithm relies on the accuracy of the edgedetection procedure it encounters
diculties with images in which regions are both noisy and have blurred or
indistinct boundaries Another interesting method within this category called
variableorder surface tting starts with a coarse segmentation of the
given image into several surfacecurvaturesign primitives for example pit
peak and ridge which are then rened by an iterative regiongrowing method
based on variableorder surface tting This method however may only be
suitable to the class of images where the image contents vary considerably
The main diculty with regionbased segmentation schemes lies in the se
lection of a homogeneity criterion Regionbased segmentation algorithms
have been proposed using statistical homogeneity criteria based on regional
feature analysis Bayesian probability modeling of images Markov
random elds and seedcontrolled homogeneity competition Seg
mentation algorithms could also rely on homogeneity criteria with respect to
gray level color texture or surface measures
Optimal thresholding
Suppose it is known a priori that the given image consists of only two principal
brightness levels with the prior probabilities P
and P
Consider the situation
where natural variations or noise modify the two gray levels to distributions
represented by Gaussian PDFs p
x and p
x where x represents the gray
level The PDF of the image gray levels is then
px P
p
x P
p
x
P
p
exp
x
P
p
exp
x
where
and
are the means of the two regions and
and
are their
standard deviations Let
Suppose that the dark regions in the image correspond to the background
and the bright regions to the objects of interest Then all pixels below a
threshold T may be considered to belong to the background and all pixels
above T may be considered as pixels belonging to the object of interest The
probability of erroneous classication is then
P
e
T P
Z
T
p
x dx P
Z
T
p
x dx
To nd the optimal threshold we may dierentiate P
e
T with respect to T
and equate the result to zero which leads to
P
p
T P
p
T
Applying this result to the Gaussian PDFs gives after taking logarithms and
some simplication the quadratic equation
AT
BT C
where
A
B
C
ln
P
P
The possibility of two solutions indicates that it may require two thresholds
to obtain the optimal threshold
If
a single threshold may be used given by
T
ln
P
P
Furthermore if the two prior probabilities are equal that is P
P
or if the
variance is zero that is the optimal threshold is equal to the average
of the two means
Thresholding using boundary characteristics The number of pixels
covered by the objects of interest to be segmented from an image is almost al
ways a small fraction of the total number of pixels in the image the graylevel
histogram of the image is then likely to be almost unimodal The histogram
may be made closer to being bimodal if only the pixels on or near the bound
aries of the object regions are considered
The selection and characterization of the edge or boundary pixels may be
achieved by using gradient and Laplacian operators as follows
bx y
if r
f
x y T
L
if r
f
x y T and r
f
x y
L
if r
f
x y T and r
f
x y
where r
f
x y is a gradient and r
f
x y is the Laplacian of the given image
fx y T is a threshold and L
L
represent three distinct gray levels
In the resulting image the pixels that are not on an edge are set to zero the
pixels on the darker sides of edges are set to L
and the pixels on the lighter
sides of edges are set to L
This information may be used not only to detect
objects and edges but also to identify the leading and trailing edges of objects
with reference to the scanning direction
See Section for a description of Otsus method of deriving the optimal
threshold for binarizing a given image see also Section for discussions
on a few other methods to derive thresholds
Regionoriented segmentation of images
Let R represent the region spanning the entire space of the given image
Segmentation may be viewed as a process that partitions R into n subregions
R
R
R
n
such that
n
i
R
i
R that is the union of all of the regions detected spans the
entire image then every pixel must belong to a region
R
i
is a connected region i n
R
i
R
j
i j i j that is the regions are disjoint
PR
i
TRUE for i n for example all pixels within a
region have the same intensity
PR
i
R
j
FALSE i j i j for example the intensities of the
pixels in dierent regions are dierent
where PR
i
is a logical predicate dened over the points in the set R
i
and
is the null set
A simple algorithm for region growing by pixel aggregation based upon the
similarity of a local property is as follows
Start with a seed pixel or a set of seed pixels
Append to each pixel in the region those of its connected or connected
neighbors that have properties gray level color etc that are similar
to those of the seed
Stop when the region cannot be grown any further
The results of an algorithm as above depend upon the procedure used to
select the seed pixels and the measures of similarity or inclusion criteria used
The results may also depend upon the method used to traverse the image
that is the sequence in which neighboring pixels are checked for inclusion
Splitting and merging of regions
Instead of using seeds to grow regions or global thresholds to separate an image
into regions one could initially consider to divide the given image arbitrarily
into a set of disjoint regions and then to split andor merge the regions using
conditions or predicates P
A general splitmerge procedure is as follows Assuming the image to
be square subdivide the entire image R successively into smaller and smaller
quadrant regions such that for any region R
i
PR
i
TRUE In other
words if PR FALSE divide the image into quadrants if P is FALSE
for any quadrant subdivide that quadrant into subquadrants Iterate the
procedure until no further changes are made or a stopping criterion is reached
The splitting technique may be represented as a quadtree Diculties could
exist in selecting an appropriate predicate P
Because the splitting procedure could result in adjacent regions that are
similar a merging step would be required which may be specied as follows
Merge two adjacent regions R
i
and R
k
if PR
i
R
k
TRUE Iterate until
no further merging is possible
Region growing using an additive tolerance
A commonly used regiongrowing scheme is pixel aggregation The
method compares the properties of spatially connected neighboring pixels with
those of the seed pixel the properties used are determined by homogeneity
criteria For intensitybased image segmentation the simplest property is
the pixel gray level The term additive tolerance level stands for the per
mitted absolute graylevel dierence between the neighboring pixels and the
seed pixel a neighboring pixel fmn is appended to the region if its abso
lute graylevel dierence with respect to the seed pixel is within the additive
tolerance level T
jfmn seedj T
Figure shows a simple example of additivetolerance region growing
using dierent seed pixels in a image The additive tolerance level used
in the example is T Observe that two dierent regions are obtained by
starting with two seeds at dierent locations as shown in Figure b and
Figure c In order to overcome this dependence of the region on the
seed pixel selected the following modied criterion could be used to determine
whether a neighboring pixel should be included in a region or not instead of
comparing the incoming pixel with the gray level of the seed the gray level of
a neighboring pixel is compared with the mean gray level called the running
mean
R
c
of the region being grown at its current stage R
c
This criterion
may be represented as
jfmn
R
c
j T
where
R
c
N
c
X X
mn R
c
fmn
where N
c
is the number of pixels in R
c
Figure d shows the result ob
tained with the runningmean algorithm by using the same additive tolerance
level as before T With the runningmean criterion no matter which
pixel is selected as the seed the same nal region is obtained in the present
example as long as the seed pixel is within the region which is the central
highlighted area in Figure d
In the simple scheme described above the seed pixel is always used to check
the incoming neighboring pixels even though most of them are not spatially
connected or close to the seed Such a regiongrowing procedure may fail when
a seed pixel is inappropriately located at a noisy pixel Shen
suggested the use of the current center pixel as the reference instead of
the seed pixel that was used to commence the regiongrowing procedure For
example the shaded area shown in Figure represents a region being
grown After the pixel C is appended to the region its connected neighbors
labeled as Ni i or connected neighbors marked as Ni i
would be checked for inclusion in the region using
jNi Cj T
The pixel C is called the current center pixel However because some of
the neighboring pixels N and N in the illustration in Figure are
FIGURE
Example of additivetolerance region growing using dierent seed pixels T
a Original image b The result of region growing shaded in black with
the seed pixel at c The result of region growing with the seed pixel at
d The result of region growing with the runningmean algorithm or
the current center pixel method using any seed pixel within the highlighted
region Figure courtesy of L Shen
already included in the region shown only N N and N in the case of
connectivity or N N N N N and N in the case of connectivity
are compared with their current center pixel C for region growing rather
than with the original seed pixel For the example shown in Figure this
procedure generates the same result as shown in Figure d independent
of the location of the seed pixel within the ROI when using the same additive
tolerance level T
FIGURE
Illustration of the concept of the current center pixel in region growing
Figure courtesy of L Shen
Region growing using a multiplicative tolerance
In addition to the sensitivity of the region to seed pixel selection with additive
tolerance region growing the additive tolerance level or absolute dierence in
gray level T is not a good criterion for region growing an additive tolerance
level of while appropriate for a seed pixel value or running mean of for
example may not be suitable when the seed pixel gray level or running
mean is at a dierent level such as or In order to address this problem
a relative dierence based upon a multiplicative tolerance level could be
employed Then the criterion for region growing could be dened as
jfmn
R
c
j
R
c
or
jfmn
R
c
j
fmn
R
c
where fmn is the gray level of the current pixel being checked for inclusion
and
R
c
could stand for the original seed pixel value the current center pixel
value or the runningmean gray level Observe that the two equations above
are comparable to the denitions of simultaneous contrast in Equations
and
The additive and multiplicative tolerance levels both determine the maxi
mum graylevel deviation allowed within a region and any deviation less than
this level is considered to be an intrinsic property of the region or to be noise
Multiplicative tolerance is meaningful when related to the SNR of a region or
image whereas additive tolerance has a direct connection with the standard
deviation of the pixels within the region or a given image
Analysis of region growing in the presence of noise
In order to analyze the performance of regiongrowing methods in the presence
of noise let us assume that the given image g may be modeled as an ideal
image f plus a noise image where f consists of a series of strictly uniform
disjoint or nonoverlapping regions R
i
i k and includes their
corresponding noise parts
i
i k Mathematically the image may
be expressed as
g f
where
f
i
R
i
i k
and
i
i
i k
A strictly uniform region R
i
is composed of a set of connected pixels fmn
at positions mn whose values equal a constant
i
that is
R
i
f mn j fmn
i
g
The set of regions R
i
i k is what we expect to obtain as the
result of segmentation Suppose that the noise parts
i
i k are
composed of white noise with zero mean and standard deviation
i
then we
have
g
i
R
i
i
i k
and
f
i
R
i
g
i
i
i k
As a special case when all the noise components have the same standard
deviation that is
k
and
k
where the symbol represents statistical similarity the image f may be de
scribed as
g
i
R
i
i k
and
f
i
R
i
g i k
Additivetolerance region growing is wellsuited for segmentation of this spe
cial type of image and an additive tolerance level solely determined by may
be used globally over the image However such special cases are rare in real
images A given image generally has to be modeled as in Equation
where multiplicativetolerance region growing may be more suitable with the
expectation that a global multiplicative tolerance level can be derived for all
of the regions in the given image Because the multiplicative tolerance level
could be made a function of
i
i
that is directly related to the SNR which
can be dened as log
i
i
dB for each individual region R
i
such a global
tolerance level can be found if
k
k
Iterative region growing with multiplicative tolerance
Simple region growing with xed additive or multiplicative tolerance may pro
vide good segmentation results with images satisfying either Equations
and or Equation However many images do not meet these con
ditions and one may have to employ dierent additive or multiplicative tol
erance levels at dierent locations in the given image This leads to the
problem of nding the appropriate tolerance level for each individual region
which the iterative multitolerance regiongrowing approach proposed by Shen
et al attempts to solve
When human observers attempt to recognize an object in an image they
are likely to use the information available in both an object of interest and
its surroundings Unlike most of the reported boundary detection methods
the HVS detects object boundaries based on not only the information around
the boundary itself such as the pixel variance and gradient around or across
the boundary but also on the characteristics of the object Shen et al
proposed using the information contained within a region as well as its
relationship with the surrounding background to determine the appropriate
tolerance level to grow a region Such information could be represented by a
set of features characterizing the region and its background With increasing
tolerance levels obtained using a certain step size it could be expected that the
values of the feature set at successive tolerance levels will be either the same
or similar which means that the corresponding regions are similar when the
regiongrowing procedure has identied an actual object or region Suppose
that the feature set mentioned above includes M features then the feature
vector V
k
at the tolerance level k may be expressed as
V
k
V
k
V
k
V
kM
T
The minimum normalized distance d
min
between the feature vectors at suc
cessive tolerance levels could be utilized to select the nal region
d
min
min
k
d
k
min
k
M
X
m
V
km
V
km
V
km
V
km
where k K and K is the number of tolerance values used
As an example the feature set could be chosen to include the coordinates of
the centroid of the region the number of pixels in the region a region shape
descriptor such as compactness and the ratio of the mean pixel value of the
foreground region to the mean pixel value of a suitably dened background
Figure shows the owchart of the iterative multitolerance regiongrow
ing algorithm The algorithm starts with a selected pixel called the seed pixel
as the rst region pixel Then the pixel value fmn of every connected
neighbor of the seed or all pixels belonging to the region at later stages of
the algorithm is checked for the condition described in Equation If
the condition is satised the pixel is included in the region This recursive
procedure is continued until no spatially connected pixel meets the condition
for inclusion in the region The outermost layer of connected pixels of the
region grown is then treated as the boundary or contour of the region
In order to permit the selection of the most appropriate tolerance level for
each region the iterative multitolerance regiongrowing procedure proposed
by Shen et al uses a range of the fractional tolerance value from
seed
to
with a step size of
seed
The specied feature vector is computed for the
region obtained at each tolerance level The normalized distance between the
feature vectors for successive tolerance levels is calculated The feature vector
with the minimum distance as dened in Equation is selected and the
corresponding region is retained as the nal object region
Figure a represents an ideal image composed of a series of strictly
uniform regions with two ROIs a relatively dark small triangle and a rel
atively bright occluded circle Two versions of the image with white noise
added are shown in Figure b SNR dB and Figure c SNR
dB As expected according to the model dened in Equation the
xed multiplicativetolerance regiongrowing method detects the two ROIs
with properly selected tolerance values for SNR dB and for
SNR dB the regions are shown highlighted in black and white in Figures
FIGURE
Flowchart of the iterative multitolerance regiongrowing algorithm Figure
courtesy of L Shen
d and e However with the composite images shown in Figure
f and g which were obtained by appending the two noisy images in parts
b and c of the gure the xedtolerance regiongrowing approach fails as
stated above The detected regions in the lower part are incorrect when the
growth tolerance level is set to be as shown in Figure f whereas
the detection of the occluded circle in the upper part fails with the toler
ance value of as shown in Figure g The iterative regiongrowing
technique automatically determines the correct growth tolerance level and
has thereby successfully detected both of the ROIs in the composite image as
shown in Figure h
Region growing based upon the human visual system
Despite its strength in automatically determining an appropriate tolerance
level for each region within an image the iterative multitolerance region
growing algorithm faces two limitations the iterative procedure is time
consuming and diculties exist in the selection of the range of the tolerance
value as well as the increment step size Shen et al used the
range of
seed
to with a step size of
seed
The motivation for dening
tolerance as being inversely related to the seed pixel value was to avoid the
possibilities of no change between regions grown at successive tolerance lev
els due to the use of too small a tolerance value and negligible increment
in terms of the general intensity of the region A possible solution to avoid
searching for the appropriate growth tolerance level is to make use of some of
the characteristics of the HVS
The multitolerance regiongrowing algorithm focuses mainly on the prop
erties of the image to be processed without consideration of the observers
characteristics The main purpose of an image processing method such as
segmentation is to obtain a certain result to be presented to human observers
for further analysis or for further computer analysis In either case the result
should be consistent with a human observers assessment Therefore eective
segmentation should be achievable by including certain basic properties of the
HVS
The HVS is a nonlinear system with a large dynamic range and a band
pass lter behavior The ltering property is characterized by the
reciprocal of the threshold contrast C
T
that is a function of both the fre
quency u and background luminance L The smallest luminance dierence
that a human observer can detect when an object of a certain size appears
with a certain background luminance level is dened as the JND which can
be quantied as
JND L C
T
A typical variation of threshold contrast as a function of the background lumi
nance known as the WeberFechner relationship is graphically depicted
in Figure The curve can be typically divided into two asymptotic re
Demonstration of the multiplicativetolerance regiongrowing algorithm
a An ideal test image b Test image with white noise SNR dB
c Test image with white noise SNR dB d Regions grown with xed
for the image in b e Regions grown with xed for the
image in c f Regions grown with xed g Regions grown with
xed h Regions grown with the iterative algorithm with adaptive
The original composite images in f h were obtained by combining
the images in b at the top and c at the bottom The detected regions are
highlighted in black for the triangle and in white for the occluded circle in
gures d h Figure courtesy of L Shen
gions the Rose de Vries region where C
T
decreases when the background
luminance increases with the relationship described as
C
T
L
and the Weber region where C
T
is independent of the background luminance
and the relationship obeys Webers law
C
T
JND
L
C
constant
FIGURE
A typical threshold contrast curve known as the WeberFechner relationship
Figure courtesy of L Shen
It is possible to determine the JND as a function of the background gray
level from psychophysical experiments In one such study conducted by
Shen a test was set up with various combinations of foreground and
background levels using an image containing a series of square boxes with
the width ranging from pixel to pixels and a xed space of pixels in
between the boxes Also included was a series of groups of four vertical
pixel lines with the width ranging from pixel to pixels and a spacing of
the same number of pixels between any two adjacent lines and a xed gap
of pixels in between the groups of lines Figure shows one such test
image with the background gray level set to be and the foreground at
Based upon the visibility or detection of up to the pixelwide square
and line group on a monitor a relationship between the JND and the back
ground gray level was obtained as shown in Figure In order to obtain
a general JND relation a large number of trials involving the participation
of a large number of subjects is necessary along with strict control of the
experimental environment Regardless Shen used the JND relationship
illustrated in Figure to develop a regiongrowing algorithm based upon
the characteristics of the HVS
FIGURE
Visual test image for determination of the JND as a function of background
gray level foreground is and background is in the displayed image
Figure courtesy of L Shen
The HVSbased regiongrowing algorithm starts with a connected neigh
borpixel grouping based upon the JND relationships of adjacent pixel gray
FIGURE
A relationship between the JND and the background gray level based upon a
psychophysical experiment Figure courtesy of L Shen
levels The JND condition is dened as
jp
p
j minfJNDp
JNDp
g
where p
and p
are two connected pixels This step is followed by the removal
of small regions dened as regions having fewer than ve pixels in the study of
Shen by merging with a connected region with the minimum mean gray
level dierence Then merging of connected regions is performed if any of two
neighboring regions meet the JND condition with p
and p
representing the
regions mean values The procedure is iterated until no neighboring region
satises the JND condition
Figure shows the results of region growing with the same test image as
in the test of the tolerancebased regiongrowing algorithms in Figure f
The HVSbased regiongrowing algorithm has successfully segmented the re
gions at the two SNR levels The method is not timeconsuming because a
JND table is used to determine the parameters required
FIGURE
Results of the HVSbased regiongrowing algorithm with the same test image
as in Figure f Figure courtesy of L Shen
Application Detection of calci cations by multitoler
ance region growing
Microcalcications are the most important and sometimes the only mam
mographic sign in early curable breast cancer Due to their sub
tlety detection and classication as benign or malignant are two major
problems Several researchers in the eld of mammographic image analy
sis have focused attention on the
detection of calcications Shen et al reported on a method to detect
and classify mammographic calcications based upon a multitolerance region
growing procedure shape analysis and neural networks The owchart of the
system is shown in Figure The calcication detection algorithm consists
of three steps
selection of seed pixels
detection of potential calcication regions and
conrmation of calcication regions
Selection of seed pixels One of the common characteristics of calci
cations in mammograms is that they are relatively bright due to the higher
Xray attenuation coecient or density of calcium as compared with other
normal breast tissues Hence a simple criterion to select seed pixels to search
for calcications could be based on the median or the mean value of the mam
mogram The scheme employed by Shen et al is as follows
Every pixel with a value greater than the median gray level of the
mammogram is identied as a potential seed pixel for the next
two steps of calcication detection The pixels identied as above
are processed in sequence by selecting the highest intensity pixel
remaining in rasterscan order as long the pixel has not been
included in any of the regions already labeled as calcications
The selection scheme is simple and eective in most cases but faces limi
tations in analyzing mammograms of dense breasts Removal of the parts of
the mammographic image outside the breast boundary see Section and
the highdensity area of the pectoral muscle see Section could be useful
preprocessing steps
Detection of potential calci cation regions Calcications could ap
pear in regions of varying density in mammograms Calcications present
within dense masses or superimposed by dense tissues in the process of ac
quisition of mammograms could present low graylevel dierences or contrast
with respect to their local background On the other hand calcications
present against a background of fat or lowdensity tissue would possess higher
dierences and contrast The iterative multitolerance regiongrowing method
presented in Section could be expected to perform well in this applica
tion by adapting to variable conditions as above
In the method proposed by Shen et al for the detection of calcica
tions the algorithm starts a regiongrowing procedure with each seed pixel
selected as above Every connected neighbor fmn of the pixels belonging
to the region is checked for the following condition
R
max
R
min
fmn R
max
R
min
FIGURE
Flowchart of a method for the detection and classication of mammographic
calcications Figure courtesy of L Shen
where R
max
and R
min
are the current maximum and minimum pixel values
of the region being grown and is the growth tolerance The fractional
tolerance value for region growing is increased from to with a step
size determined as the inverse of the seedpixels gray level A feature vector
including compactness dened as c
area
perimeter
see Section for
details the x y coordinates of the centroid and the size or area in number
of pixels is calculated for the region obtained at each tolerance level The
normalized distance between the feature vectors for successive tolerance levels
is computed as given by Equation The feature set with the minimum
distance is selected as the nal set and the corresponding region considered
to be a potential calcication region
Con rmation of calci cation regions Each potential calcication re
gion detected is treated as a calcication region only if the size S in pixels
and contrast C computed as in Equation of the region at the nal level
meet the following conditions
S
and
C
The upper limit on the area corresponds to about mm
with a pixel
resolution of m The background region required to compute C is formed
by using pixels circumscribing the region contour to a thickness of pixels
The contrast threshold of was selected based upon another study on seg
mentation and analysis of calcications
Examples Two examples of the detection of calcications by the method
described above are presented in Figures and Figure a and
Figure a show two sections of mammograms of size pixels
with benign calcications and malignant calcications respectively Figure
b and Figure b show the same mammogram sections with the
contours of the calcication regions extracted by the algorithm as described
above
Sections of size and pixels
of four typical mammograms from complete images of up to
pixels with biopsyproven calcications were used in the study by Shen et
al Two of the sections had a total of benign calcications whereas
the other two contained malignant calcications Based upon visual
inspection by a radiologist the detection rates of the multitolerance region
growing algorithm were ! with false calcications and ! with
false calcications for the benign and malignant mammograms respectively
Bankman et al compared their hillclimbing segmentation algorithm
for detecting microcalcications with the multitolerance regiongrowing algo
rithm described above Their results showed that the two algorithms have sim
ilar discrimination powers based on an ROC analysis conducted with six mam
mograms containing clusters with a total of calcications Bankman
et al stated that The multitolerance algorithm provides a good solution to
avoid the use of statistical models local statistic estimators and the manual
selection of thresholds However the cost is multiple segmentations of the
same structure and computation of features during the segmentation of each
structure The segmented regions were comparable in many cases
Details on further analysis of the calcications detected as above are pro
vided in Sections and See Section for another method for the
detection of calcications
Application Detection of calci cations by linear pre
diction error
The simple seed selection method used by Shen et al and described
in Section encounters limitations in the case of calcications present
in or superimposed by dense breast tissue Serrano et al proposed a
method to detect seed pixels for region growing based upon the error of linear
prediction The D linear prediction error computation method proposed by
Kuduvalli and Rangayyan was used to compute the prediction
error directly from the image data without the need to compute the prediction
coecients see Section for details
The method proposed by Serrano et al starts with a preltering step
Considering the small size of microcalcications a lowpass lter with a wide
kernel could be expected to remove them from the image while conserving a
background of high density in the image Conversely a highpass lter may be
used to detect microcalcications Serrano et al employed a highpass lter
specied as hmn gmn where gmn was a lowpass Gaussian
function with a variance of pixels a lter kernel size of pixels
was chosen In the next step a D linear prediction error lter
was applied to the output of the highpass lter A pixel was selected
as a seed for the multitolerance regiongrowing algorithm if its prediction
error was greater than an experimentally determined threshold This was
based upon the observation that a microcalcication can be seen as a point of
nonstationarity in an approximately homogeneous region or neighborhood in
a mammogram such a pixel cannot be predicted well by the linear predictor
and hence leads to a high error
The detection algorithm was tested with three mammograms containing a
total of microcalcications of dierent nature and diagnosis The results
obtained with the algorithm were examined by a radiologist who determined
the accuracy of the detection Figure shows a segment of a mammogram
with several calcications the seed pixels identied by thresholding the pre
diction error and the regions obtained by application of the multitolerance
regiongrowing algorithm to the seed pixels In comparison with the algorithm
of Shen et al the detection accuracy of the method of Serrano et al
was higher with a smaller number of false detections for the images tested
a
b
FIGURE
Mammogram section with benign calcications a Original image b Im
age with the contours of the calcication regions detected The section shown
is of size pixels approximately cm cm out of the
full matrix of pixels of the complete mammogram Repro
duced with permission from L Shen RM Rangayyan and JEL Desautels
Detection and classication of mammographic calcications International
Journal of Pattern Recognition and Articial Intelligence
c
World Scientic Publishing Co
a
b
FIGURE
Mammogram section with malignant calcications a Original image
b Image with the contours of the calcication regions detected The sec
tion shown is of size pixels approximately cm cm out
of the full matrix of pixels of the complete mammogram Repro
duced with permission from L Shen RM Rangayyan and JEL Desautels
Detection and classication of mammographic calcications International
Journal of Pattern Recognition and Articial Intelligence
c
World Scientic Publishing Co
although the detection capability was diminished that is more calcications
were missed by the predictionerror method
FIGURE
a Mammogram section with malignant calcications pixels with
a resolution of m b Seed pixels detected by thresholding the predic
tion error marked in black c Image with the contours of the calcication
regions detected by region growing from the seed pixels in b Reproduced
with permission from C Serrano JD Trujillo B Acha and RM Ran
gayyan Use of D linear prediction error to detect microcalcications in
mammograms CDROM Proceedings of the II Latin American Congress on
Biomedical Engineering Havana Cuba May
c
Cuban Society
of Bioengineering
Fuzzysetbased Region Growing to Detect Breast
Tumors
Although mammography is being used for breast cancer screening the anal
ysis of masses and tumors on mammograms is at times dicult because
developing signs of cancer may be minimal or masked by superimposed tis
sues Additional diagnostic procedures may be recommended when the origi
nal mammogram is equivocal
Computeraided image analysis techniques have the potential to improve
the diagnostic accuracy of mammography and reduce the use of adjunctive
procedures and morbidity as well as healthcare costs Computer analysis can
facilitate the enhancement detection characterization and quantication of
diagnostic features such as the shapes of calcications and masses growth
of tumors into surrounding tissues and the distortion caused by developing
densities The annotation of mammograms with objective measures may assist
radiologists in diagnosis
Computeraided detection of breast masses is a challenging problem requir
ing sophisticated techniques due to the low contrast and poor denition of
their boundaries Classical segmentation techniques attempt to dene pre
cisely the ROI such as a calcication or a mass Shen et al proposed
thresholding and multitolerance region growing methods for the detection of
potential calcication regions and extraction of their contours see Sections
and Karssemeijer Laine et al and Miller and Ram
sey proposed methods for tumor detection based on scalespace analysis
Zhang et al proposed an automated detection method for the initial iden
tication of spiculated lesions based on an analysis of mammographic texture
patterns Matsubara et al described an algorithm based on an adaptive
thresholding technique for mass detection Kupinski and Giger presented
two methods for segmenting lesions in digital mammograms a radialgradient
indexbased algorithm that considers both the graylevel information and a
geometric constraint and a probabilistic approach However dening criteria
to realize precisely the boundaries of masses in mammograms is dicult The
problem is compounded by the fact that most malignant tumors possess fuzzy
boundaries with a slow and extended transition from a dense core region to
the surrounding tissues For detailed reviews on the detection and analy
sis of breast masses refer to Rangayyan et al and Mudigonda et
al See Sections and for related discussions
An alternative to address the problem of detecting breast masses is to rep
resent tumor or mass regions by fuzzy sets The most popular algorithm
that uses the fuzzyset approach is the fuzzy Cmeans algorithm
The fuzzy Cmeans algorithm uses iterative optimization of an objective
function based on weighted similarity measures between the pixels in the im
age and each cluster center The segmentation method of Chen and Lee
uses fuzzy Cmeans as a preprocessing step in a Bayesian learning paradigm
realized via the expectationmaximization algorithm for edge detection and
segmentation of calcications and masses in mammograms However their
nal result is based on classical segmentation to produce crisp boundaries
Sameti and Ward proposed a lesion segmentation algorithm using fuzzy
sets to partition a given mammogram Their method divides a mammogram
into two crisp regions according to a fuzzy membership function and an it
erative optimization procedure to minimize an objective function If more
than two regions are required the algorithm can be applied to each region
obtained in the preceding step using the same procedure The authors pre
sented results of application of the method to mammograms with four levels
of segmentation
Guliato et al proposed two segmentation methods that incorporate
fuzzy concepts The rst method determines the boundary of a tumor or mass
by region growing after a preprocessing step based on fuzzy sets to enhance the
ROI The second segmentation method is a fuzzy regiongrowing method that
takes into account the uncertainty present around the boundaries of tumors
These methods are described in the following sections
Preprocessing based upon fuzzy sets
A mass or tumor typically appears on a mammogram as a relatively dense
region whose properties could be characterized using local density gradient
texture and other measures A set of such local properties could be used to
dene a feature vector of a mass ROI andor a pixel belonging to the ROI
Given a feature vector a pixel whose properties are similar to those repre
sented by the feature vector of the mass could be assigned a high intensity
If the properties do not match the pixel intensity could be made low At
the end of such a process the pixels in and around the ROI will be displayed
according to their degree of similarity with respect to the features of the mass
ROI
A fuzzy set may be dened by assigning to each element considered from
the universal set " a value representing its grade of membership in the fuzzy
set The grade corresponds to the degree with which the element
is similar to or compatible with the concept represented by the fuzzy set Let
# " L be a membership function that maps " into L where L denotes
any set that is at least partially ordered The most commonly used range of
values for membership functions is the unit real interval Crisp sets can
be seen as a particular case of fuzzy sets where # " f g that is the
range includes only the discrete values and
The enhancement and subsequent detection of an ROI may be achieved
by dening an appropriate membership function that evaluates the similarity
between the properties of the pixel being considered and those of the ROI
itself given by the feature vector In this procedure the original image is
mapped to a fuzzy set according to the membership function which
assigns a membership degree equal to to those pixels that possess the
same properties as the mass ROI
represents the degree of similarity between the features of the mass ROI
and those of the pixel being considered
exhibits symmetry with respect to the dierence between the features
of the ROI and those of the pixel being considered and
decreases monotonically from to
Guliato et al considered the mean intensity of a seed region identied
by the user as the ROI feature A membership function with the character
istics cited above illustrated in Figure is given by the function
#p
j AB j
where p is the pixel being processed A is the feature vector of the mass gray
level B is the feature vector of the pixel being analyzed and denes the
opening of the membership function For large the opening is narrow and
the functions behavior is strict for small the opening is wide and the
function presents a more permissive behavior
FIGURE
Fuzzy membership function for preprocessing Reproduced with permission
from D Guliato RM Rangayyan WA Carnielli JA Zuo and JEL
Desautels Segmentation of breast tumors in mammograms using fuzzy sets
Journal of Electronic Imaging
c
SPIE and IS$T
The fuzzy set obtained by the method described above represents pixels
whose properties are close to those of the mass with a high membership degree
the opposite case results in a low membership degree The membership degree
may be used as a scale factor to obtain gray levels and display the result as
an image The contrast of the ROI in the resulting image depends upon the
parameter
Figures a and b show a pixel portion of a mammogram
with a spiculated malignant tumor and the result of fuzzysetbased prepro
cessing with respectively It is seen from the image in Figure
b that the pixels in the tumor region the bright area in the upperleft
part of the image have higher values than the pixels in other parts of the
image indicating a higher degree of similarity with respect to the ROI or
seed region The membership values decrease gradually across the boundary
of the tumor as expected due to the malignant nature of the tumor in this
particular case Note however that a few other spatially disconnected regions
on the righthand side of the image also have high values these regions can
be eliminated by further processing as described in Section
Figure a
Fuzzy segmentation based upon region growing
Region growing is an image segmentation technique that groups pixels or
subregions into larger regions according to a similarity criterion Statistical
measures provide good tools for dening homogeneous regions The success of
image segmentation is directly associated with the choice of the measures and
a suitable threshold In particular mean and standard deviation measures are
often used as parameters to control region growing however these measures
are in uenced by extreme pixel values As a consequence the nal shape
of the region grown depends upon the strategy used to traverse the image
Figure b
Figure c
d
FIGURE
a A pixel portion of a mammogram with a spiculated malignant
tumor Pixel size m b Fuzzysetbased ROI enhancement with
c Contour extracted white line by region growing with the result in
b The black line represents the boundary drawn by a radiologist shown for
comparison threshold d Result of fuzzy region growing
with the image in a with
max
CV
max
The
contour drawn by the radiologist is superimposed for comparison Reproduced
with permission from D Guliato RM Rangayyan WA Carnielli JA Zuo
and JEL Desautels Segmentation of breast tumors in mammograms using
fuzzy sets Journal of Electronic Imaging
c
SPIE
and IS$T
Furthermore the algorithm could present unstable behavior for example dif
ferent pixels with the same values that are rejected at an earlier stage may be
accepted later on in the regiongrowing method It is also possible that the
stopping condition is not reached when the gray level in the image increases
slowly in the same direction as that of the traversal strategy Besides tradi
tional regiongrowing methods represent the ROI by a classical set dening
precisely the regions boundary In such a case the transition information is
lost and the segmentation task becomes a critical stage in the image analysis
system In order to address these concerns Guliato et al presented
two image segmentation methods the rst based on classical region growing
with the fuzzyset preprocessed image described in the following paragraphs
and the second based on fuzzy region growing using statistical measures in
homogeneity criteria described in Section
The pixel values in the fuzzyset preprocessed image represent the member
ship degrees of pixels with respect to the ROI as dened by the seed region To
perform contour extraction the regiongrowing algorithm needs a threshold
value and a seed region that lies inside the ROI The regiongrowing process
starts with the seed region Fourconnected neighboring pixels that are above
the threshold are labeled as zero the neighbors of the pixels labeled as zero
are inspected and the procedure continued If the connected pixel is less than
the threshold it is labeled as one indicating a contour pixel and its neigh
borhood is not processed The recursive process continues until all spatially
connected pixels fail the test for inclusion in the region A postprocessing step
is included to remove isolated pixels and regions that lie within the outermost
contour
The algorithm is simple and easy to implement and will always produce
closed contours The method was evaluated with a number of synthetic test
images as well as medical images such as CT and nuclear medicine images
and produced good results even in the presence of high levels of noise
Examples Figure shows the results of the method with a synthetic
image for three representative combinations of parameters The three results
exhibit a good degree of similarity and illustrate the robustness of the method
in the presence of noise
Figure c shows the contour extracted for the mammogram in part
a of the same gure Figure a shows a part of a mammogram with
a circumscribed benign mass part b of the gure shows the corresponding
enhanced image Figure c shows the contour obtained the image
is superimposed with the contour obtained by region growing in white the
contour in black is the boundary drawn independently by an experienced
radiologist shown for comparison
Results of application to mammograms Guliato et al tested
their method with mammograms including malignant tumors and
benign masses and observed good agreement between the contours given by
the method and those drawn independently by a radiologist The seed re
gion and threshold value were selected manually for each case the threshold
a b
c d
FIGURE
Illustration of the eects of seed pixel and threshold selection on fuzzyset
preprocessing and region growing a Original image pixels with
additive Gaussian noise with and SNR Results with b seed
pixel and threshold c seed pixel and threshold
d seed pixel and threshold Reproduced with permission
from D Guliato RM Rangayyan WA Carnielli JA Zuo and JEL
Desautels Segmentation of breast tumors in mammograms using fuzzy sets
Journal of Electronic Imaging
c
SPIE and IS$T
Figure a
values varied between and for the images used The same value of
the membership function parameter was used to process all of the
images in the study It was observed that the result of segmentation depended
upon the choice of the seed to start region growing and the threshold Au
tomatic selection of the seed pixel or region and the threshold is a dicult
problem that was not addressed in the study It was observed that the thresh
old could possibly be derived as a function of the statistics such as the mean
and standard deviation of the fuzzyset preprocessed image
Measure of fuzziness In order to compare the results obtained by seg
mentation with the contours of the masses drawn by the radiologist Guliato
et al developed a method to aggregate the segmented region with
the reference contour The procedure can aggregate not only two contours
but also a contour with a fuzzy region and hence is more general than clas
sical intersection The method uses a fuzzy fusion operator that generalizes
classical intersection of sets producing a fuzzy set that represents the agree
ment present among the two inputs see Section for details The result
of fusion was evaluated by a measure of fuzziness computed as
fX
P
pX
j #p j
j X j
Figure b
Figure c
d
FIGURE
a A pixel portion of a mammogram with a circumscribed be
nign mass Pixel size m b Fuzzysetbased ROI enhancement with
c Contour extracted white line by region growing with the
result in b The black line represents the boundary drawn by a radiologist
shown for comparison threshold d Result of fuzzy
region growing with the image in a with
max
CV
max
The contour drawn by the radiologist is superimposed for compar
ison Reproduced with permission from D Guliato RM Rangayyan WA
Carnielli JA Zuo and JEL Desautels Segmentation of breast tumors
in mammograms using fuzzy sets Journal of Electronic Imaging
c
SPIE and IS$T
where X is the result of aggregation and #p is the degree of membership of
the pixel p The denominator in the expression above normalizes the measure
with respect to the area of the result of fusion resulting in a value in the
range with zero representing perfect agreement and unity indicating no
intersection between the two inputs
The values of the measure of fuzziness obtained for the mammograms in
the study were in the range with the mean and standard deviation
being and respectively The measure of fuzziness was less than
for out of the cases In most cases where the measure of fuzziness was
greater than the segmented region was smaller than but contained within
the region indicated by the contour drawn by the radiologist Regardless of
the agreement in terms of the measure of fuzziness it was argued that for
a spiculated lesion there is no denite number of spicules that characterizes
the lesion as malignant The method captured the majority of the spicules in
the cases analyzed providing sucient information for diagnosis according
to the analysis of the results performed by an expert radiologist
Assessment of the results by pattern classi cation In order to de
rive a parameter for discriminating between benign masses and malignant
tumors the following procedure was applied by Guliato et al A
morphological erosion procedure with a square structuring element of size
equal to ! of the shorter dimension of the smallest rectangle containing the
contour was applied to the contour so that the core of the ROI was separated
from the boundary A parameter labeled as DCV was computed from the
fuzzyset preprocessed image by taking the dierence between the coecient
of variation CV of the entire ROI and that of the core of the ROI A high
value of DCV represents an inhomogeneous ROI which could be indicative
of a malignant tumor The probability of malignancy based upon DCV was
computed using the logistic regression method see Section for details
the result is illustrated in Figure Several cut points were analyzed with
the curve the cut point of resulted in all benign masses and out of
the malignant tumors being correctly classied yielding a high specicity
of but a low sensitivity of
Fuzzy region growing
Guliato et al also proposed a fuzzy regiongrowing algorithm to obtain
mass regions in mammograms In this method an adaptive similarity crite
rion is used for region growing with the mean and the standard deviation
of the pixels in the region being grown as control parameters The region
is represented by a fuzzy set to preserve the transition information around
boundary regions
The algorithm starts with a seed region that lies inside the ROI and spreads
by adding to the region connected pixels that have similar properties The
homogeneity of the region is evaluated by calculating the mean standard
deviation and the coecient of variation CV
FIGURE
The probability of malignancy vertical axis derived from the parameter
DCV horizontal axis Reproduced with permission from D Guliato RM
Rangayyan WA Carnielli JA Zuo and JEL Desautels Segmentation
of breast tumors in mammograms using fuzzy sets Journal of Electronic
Imaging
c
SPIE and IS$T
Let
max
CV
max
and be the control parameters for region growing
max
species the maximum allowed dierence between the value of the
pixel being analyzed and the mean of the subregion already grown CV
max
indicates the desired degree of homogeneity between two subregions denes
the opening of the membership function Let p be the next pixel to be analyzed
and Ip be the value of p The segmentation algorithm is executed in two
steps
j Ip j
max
If this condition is not satised then the pixel is
labeled as rejected If the condition is satised p is temporarily added
to the subregion and
new
and
new
are calculated
j
new
new
j CV
max
If the condition is satised then p must de
nitely be added to the subregion and labeled as accepted and and
must be updated that is
new
and
new
If the condition is
not satised p is added to the subregion with the label accepted with
restriction and and are not modied
The second step given above analyzes the distortion that the pixel p can
produce if added to the subregion At the beginning of the process the region
includes all the pixels in the seed region and the standard deviation is set
to zero While the standard deviation of the region being grown is zero a
specic procedure is executed in the second step j
new
new
j CV
max
The parameter CV
max
works as a lter that avoids the possibility that the
mean and standard deviation measures suer undesirable modication during
the regiongrowing process Furthermore the algorithm processes pixels in
expanding concentric squares around the seed region evaluating each pixel
only once These steps provide stability to the algorithm
The membership function that maps the pixel values of the region resulting
from the preceding procedure to the unit interval could be based upon
the mean of the region Pixels that are close to the mean will have a high
membership degree and in the opposite case a low membership degree The
desirable characteristics of the membership function are
the membership degree of the seed pixel or region must be
the membership degree of a pixel labeled as rejected must be
the membership function must be as independent of the seed pixel or
region as possible
the membership degree must represent the proximity between a pixel
labeled as accepted or accepted with restriction and the mean of the
resulting region
the function must be symmetric with respect to the dierence between
the mean and the pixel value and
the function must decrease monotonically from to
The membership function # used by Guliato et al is illustrated in
Figure where a j mean seed region j and b
max
The value
of a pixel p is mapped to the fuzzy membership degree #p as follows
if j Ip j a then #p
else if j Ip j b then #p
else #p
jIp j
The method was tested on several synthetic images with various levels of
noise Figure illustrates three representative results of the method with
a synthetic image and dierent seed pixels The results do not dier signi
cantly indicating the low eect of noise on the method
Results of application to mammograms The fuzzy region for the
malignant tumor shown in Figure a is illustrated in part d of the
same gure Figure d shows the fuzzy region obtained for the benign
mass shown in part a of the same gure
An interactive graphical interface was developed by Guliato et al
using an objectoriented architecture with controller classes Some of the fea
tures of the interface are fast and easy upgradability portability and threads
FIGURE
Fuzzy membership function for region growing where a j mean seed region
j and b
max
Reproduced with permission from D Guliato RM
Rangayyan WA Carnielli JA Zuo and JEL Desautels Segmentation
of breast tumors in mammograms using fuzzy sets Journal of Electronic
Imaging
c
SPIE and IS$T
to support parallelism between tasks The interface integrates procedures to
detect contours using fuzzy preprocessing and region growing extract fuzzy
regions using fuzzy region growing compute statistical parameters and clas
sify masses and tumors as benign or malignant The interface also provides
access to basic image processing procedures including zooming in or out l
ters histogram operations the Bezier method to manipulate contours and
image format conversion
Guliato et al applied the fuzzy regiongrowing method to test
images maintaining the same values of and CV
max
and
varying only the parameter
max
The values of the parameters were se
lected by comparing the results of segmentation with the contours drawn by
a radiologist The
max
parameter ranged from to for the masses
and tumors analyzed
The fuzzy regions obtained for the mammograms were compared objec
tively with the corresponding contours drawn by the radiologist by computing
the measure of fuzziness as in Equation The values were distributed over
the range with the mean and standard deviation being and
respectively The measure of fuzziness was smaller than in of
the cases analyzed Regardless of this measure of agreement it was found
that the fuzzy regions segmented contained adequate information to facilitate
discrimination between benign masses and malignant tumors as described
next
a b
c d
FIGURE
Illustration of the eects of seed pixel selection on fuzzy region growing
a Original image pixels with Gaussian noise with and
SNR Results with b seed pixel
max
CV
max
c seed pixel
max
CV
max
d seed pixel
max
CV
max
Re
produced with permission from D Guliato RM Rangayyan WA Carnielli
JA Zuo and JEL Desautels Segmentation of breast tumors in mam
mograms using fuzzy sets Journal of Electronic Imaging
c
SPIE and IS$T
Assessment of the results by pattern classi cation In order to de
rive parameters for pattern classication Guliato et al analyzed the
characteristics of a fuzzy ribbon dened as the connected region whose pixels
possess membership degrees less than unity and separate the tumor core from
the background as illustrated in Figure Shape factors of mass contours
as well as measures of edge sharpness and texture have been proposed for the
purpose of classication of breast masses see Sections
and for related discussions However important infor
mation is lost in analysis based on crisply dened contours the uncertainty
present in andor around the ROI is not considered Guliato et al evaluated
the potential use of statistical measures of each segmented fuzzy region and
of its fuzzy ribbon as tools to classify masses as benign or malignant Observe
that the fuzzy ribbon of the malignant tumor in Figure a contains more
pixels with low values than that of the benign mass in part b of the same
gure This is due to the fact that in general malignant tumors possess
illdened boundaries whereas benign masses are wellcircumscribed Based
upon this observation Guliato et al computed the coecient of variation
CV
fr
of the membership values of the pixels lying only within the fuzzy rib
bon and the ratio
fr
of the number of pixels with membership degree less
than to the total number of pixels within the fuzzy ribbon It was expected
that the fuzzy ribbons of malignant tumors would possess higher CV
fr
and
fr
than those of benign masses
In pattern classication experiments discrimination between benign masses
and malignant tumors with the parameter
fr
had no statistical signicance
The probability of malignancy curve based upon CV
fr
computed using the
logistic regression method see Section for details is illustrated in Fig
ure The cut point of resulted in the correct classication of out
of malignant tumors and out of benign masses processed leading to
a sensitivity of and a specicity of
The fuzzy segmentation techniques described above represent the ROI by
fuzzy sets instead of crisp sets as in classical segmentation The results of
the fuzzy approach agree well with visual perception especially at transitions
around boundaries The methods allow the postponement of the crisp decision
to a higher level of image analysis For further theoretical notions related to
fuzzy segmentation and illustrations of application to medical images see
Udupa and Samarasekera and Saha et al
Detection of Objects of Known Geometry
Occasionally images contain objects that may be represented in an analyti
cal form such as straightline segments circles ellipses and parabolas For
Figure a
example the edge of the pectoral muscle appears as an almoststraight line
in MLO mammograms benign calcications and masses appear as almost
circular or oval objects in mammograms several types of cells in pathology
specimens have circular or elliptic boundaries and some may have nearly rect
angular shapes and parts of the boundaries of malignant breast tumors may
be represented using parabolas The detection modeling and characteriza
tion of objects as above may be facilitated by prior knowledge of their shapes
In this section we shall explore methods to detect straight lines and circles
using the Hough transform
The Hough transform
Hough proposed a method to detect straight lines in images based upon
the representation of straight lines in the image x y space using the slope
intercept equation
y m x c
where m is the slope and c is the position where the line intercepts the y axis
see Figure In the Hough domain or space straight lines are character
ized by the pair of parameters m c the Hough space is also known as the
parameter space A disadvantage of this representation is that both m and
c have unbounded ranges which creates practical diculties in the computa
b
FIGURE
The fuzzy ribbons of a the malignant tumor in Figure a and b the
benign mass in Figure a Reproduced with permission from D Gu
liato RM Rangayyan WA Carnielli JA Zuo and JEL Desautels
Segmentation of breast tumors in mammograms using fuzzy sets Journal
of Electronic Imaging
c
SPIE and IS$T
FIGURE
The probability of malignancy vertical axis derived from the parameter CV
fr
horizontal axis Reproduced with permission from D Guliato RM Ran
gayyan WA Carnielli JA Zuo and JEL Desautels Segmentation of
breast tumors in mammograms using fuzzy sets Journal of Electronic Imag
ing
c
SPIE and IS$T
tional representation of the m c space In order to overcome this limitation
Duda and Hart proposed the representation of straight lines using the
normal parameters as
x cos y sin
see Figure This representation has the advantage that is limited to the
range or and is limited by the size of the given image The
origin may be chosen to be at the center of the given image or at any other
convenient point the limits of the parameters are aected by the choice
of the origin
A procedure for the detection of straight lines using this parametric rep
resentation is described in the next section The Hough transform may be
extended for the detection of any curve that may be represented in a para
metric form see Rangayyan and Krishnan for an application of the
Hough transform for the identication of linear sinusoidal and hyperbolic
frequencymodulated components of signals in the timefrequency plane
Detection of straight lines
Suppose we are given a digital image that contains a straight line Let the
pixels along the line be represented as fxn yng n N
where N is the number of pixels along the line It is assumed that the image
has been binarized such that the pixels that belong to the line have the value
FIGURE
Parametric representation of a straight line in three coordinate systems x y
m c and
and all other pixels have the value It is advantageous if the line is
onepixel thick otherwise several lines could exist within a thick line
If the normal parameters of the line are
all pixels along the line
satisfy the relationship
xn cos
yn sin
For a given pixel fxn yng this represents a sinusoidal curve in the
parameter space it follows that the curves for all the N pixels intersect at
the point
The following properties of the above representation follow
A point in the x y space corresponds to a sinusoidal curve in the
parameter space
A point in the space corresponds to a straight line in the x y
space
Points lying on the same straight line in the x y space correspond to
curves through a common point in the parameter space
Points lying on the same curve in the parameter space correspond to
lines through a common point in the x y space
Based upon the discussion above a procedure to detect straight lines is as
follows
Discretize the parameter space into bins by quantizing and
as
k
k K and
l
l L the bins
are commonly referred to as accumulator cells Suitable limits may be
imposed on the ranges of the parameters
For each point in the given image that has a value of increment by
each accumulator cell in the space that satises the relationship
xn cos yn sin Note that exact equality needs to be trans
lated to a range of acceptance depending upon the discretization step
size of the parameter space
The coordinates of the point of intersection of all the curves in the
parameter space provide the parameters of the line This point will
have the highest count in the parameter space
The procedure given above assumes the existence of a single straight line in
the image If several lines exist there will be the need to search for all possible
points of intersection of several curves or the local maxima Note that the
count in a given accumulator cell represents the number of pixels that lie on
a straight line or several straightline segments that have the corresponding
parameters A threshold may be applied to detect only lines that have a
certain minimum length number of pixels All cells in the parameter space
that have counts above the threshold may be taken to represent straight lines
or segments with the corresponding values and numbers of pixels
Examples Figure a shows an image with a single straight line
represented by the parameters
o
The limits of the x and y
axes are with the origin at the center of the image The Hough transform
of the image is shown in part b of the gure in the parameter space
The maximum value in the parameter space occurs at
o
An image containing two lines with
o
and
o
is
shown in Figure a along with its Hough transform in part b of the
gure The value of was considered to be negative for normals to lines
extending below the horizontal axis x in the image with the origin at
the center of the image the range of was dened to be
o
It is also
possible to maintain to be positive with the range of extended to
o
The parameter space clearly demonstrates the expected sinusoidal patterns
as well as two peaks at the locations corresponding to the parameters of the
two lines present in the image Observe that the intensity of the point at
the intersection of the sinusoidal curves for the second line the lower of the
two bright points in the parameter space is less than that for the rst line
re ecting its shorter length
The application of the Hough transform to detect the pectoral muscle in
mammograms is described in Section See Section for further discus
sion on the Hough transform and for a modication of the Hough transform
by inclusion of the Radon transform
a b
FIGURE
a Image with a straight line with
o
The limits of the x and y
axes are with the origin at the center of the image b Hough transform
parameter space for the image The display intensity is log accumulator
cell value The horizontal axis represents
o
the vertical axis
represents
Detection of circles
For example all points along the perimeter of a circle of radius c centered at
x y a b satisfy the relationship
x a
y b
c
Any circle is represented by a single point in the D a b c parameter space
The points along the perimeter of a circle in the x y plane describe a right
circular cone in the a b c parameter space The algorithm for the detection
of straight lines described in Section may be easily extended for the
detection of circles using this representation
Example A circle of radius pixels centered at x y in a
image is shown in Figure The Hough parameter a b c space
a b
FIGURE
a Image with two straight lines with
o
and
o
The
limits of the x and y axes are with the origin at the center of the image
b Hough transform parameter space for the image The display intensity is
log accumulator cell value The horizontal axis represents
o
the vertical axis represents
of the circle is shown in Figure The parameter space demonstrates a
clear peak at a b c as expected
FIGURE
A image with a circle of radius pixels centered at x y
An image derived from the scartissue collagen ber image in Figure b
is shown in Figure a The image was binarized with a threshold equal
to of its maximum value In order to detect the edges of the objects in
the image an image was created with the value zero at all pixels having the
value zero and also having all of their connected pixels equal to zero in the
binarized image The same step was applied to all pixels with the value of one
All remaining pixels were assigned the value of one Figure b shows
the result that depicts only the edges of the bers which are nearly circular
in crosssection Observe that some of the object boundaries are incomplete
and not exactly circular The Hough parameter a b c space of the image is
shown in Figure for circles of radius in the range pixels Several
distinct peaks are visible in the planes for radius values of and pixels
the locations of the peaks give the coordinates of the centers of the circles that
are present in the image and their radii The parameter space could be further
processed and thresholded to detect automatically the peaks present which
relate to the circles in the image Prior knowledge of the range of possible
radius values could assist in the process
Figure c shows the parameter space plane for a radius of pixels
superimposed on a reversed version of the edge image in Figure b the
edges are shown in black The composite image demonstrates clearly that
the peaks in the parameter space coincide with the centers of nearly circular
objects with an estimated radius of pixels no peaks or high values are
present within smaller or larger objects A similar composite image is shown
in part d of the gure for a radius of pixels Similar results are shown in
Figures and for another TEM image of a normal ligament sample
FIGURE
Hough parameter a b c space of the circle image in Figure Each image
is of size pixels and represents the range of the a b coordinates of
the center of a potential circle The series of images represents the various
planes of the a b c parameter space with c going left to right
and top to bottom representing the radius of potential circles The intensity
of the parameter space values has been enhanced with the log operation The
maximum value in the Hough space is located at the center of the plane for
c
Observe that the Hough transform works well even when the given image has
incomplete and slightly distorted versions of the pattern being represented
Frank et al performed an analysis of the diameter distribution of col
lagen bers in rabbit ligaments Scartissue samples from injured and healing
ligaments were observed to have an almostuniform distribution of ber di
ameter in the range nm whereas normal samples were observed to
have a wider range of diameter with an average value of about nm
Methods for the Improvement of Contour or Region
Estimates
It is often the case that the contour of an ROI provided by an image processing
technique does not satisfy the user The user may wish to impose conditions
such as smoothness on the contour It may also be desirable to have the
result authenticated by an expert in the eld of application In such cases
the need may arise to modify the contour A few techniques that may assist
in this task are summarized below
Polygonal and parabolic models The contour on hand may be seg
mented into signicant parts by selecting control points or by detecting its
points of in ection see Section The contour in its entirety or its
segments may then be approximated by parametric curves such as a poly
gon see Section or parabolas see Section Minor arti
facts details or errors in the contour are removed by the parametric models
to the extent permitted by the specied tolerance or error
Bspline and Bezier curves Bezier polynomials may be used to t
smooth curves between guiding points that are specied in an interactive man
ner This approach may be used to modify a part of a contour without
aecting the rest of the contour In the case where several guiding or critical
points are available around an ROI and a smooth curve is desired to connect
them Bsplines may be used to derive the interpolating points The
guiding points may also be used for compact representation of the contour
Active contour models or snakes Kass et al proposed active
contour models or snakes that allow an initial contour to reshape and mold
itself to a desired ROI based upon constraints related to the derivatives of
the contour image gradient and energy functionals A snake is an energy
minimizing spline that is guided by external constraint forces in uenced by
imagebased forces that pull it toward features such as edges and constrained
by internal spline forces that impose piecewise smoothness The initial contour
may be provided by an image processing technique or could be provided by
the user in the form of a simple polygonal circular or oval contour within
the ROI The active contour model determines the true boundary of the
a b
c d
FIGURE
a TEM image showing collagen bers in crosssection a part of the image
in Figure b The image is of size pixels b Edges extracted
from the image in a c Negative version of the image in b overlaid with
times the c plane of the Hough transform parameter space d Same
as in c but with the c plane See also Figure
FIGURE
Hough parameter a b c space of the image in Figure b Each image
is of size pixels and represents the range of the a b coordinates of
the center of a potential circle The series of images represents the various
planes of the a b c parameter space with c going left to right
and top to bottom representing the radius of potential circles The intensity
of the parameter space values has been enhanced with the log operation
a b
c b
FIGURE
a TEM image showing collagen bers in crosssection a part of the image
in Figure a The image is of size pixels b Edges extracted
from the image in a c Negative version of the image in b overlaid with
times the c plane of the Hough transform parameter space d Same
as in c but with the c plane See also Figure
FIGURE
Hough parameter a b c space of the image in Figure b Each image
is of size pixels and represents the range of the a b coordinates
of the center of a potential circle The series of images represents the various
planes of the a b c parameter space with c going left to right
and top to bottom representing the radius of potential circles The intensity
of the parameter space values has been enhanced with the log operation
ROI that satises the conditions imposed The use of active contour models
to obtain the breast boundary in mammograms is described in Section
The live wire Falc%ao et al argued that there will continue
to be situations where automatic image segmentation methods fail and pro
posed a usersteered segmentation paradigm labeled as the live wire Their
guiding principles were to provide eective control to the user over the seg
mentation process during execution and to minimize the users time required
in the process of segmentation In the livewire approach the user initially
species a point on the boundary of the ROI using a cursor Then as the user
moves the cursor an optimal path connecting the initial point to the current
cursor position is computed and displayed in real time Optimal paths are
determined by computing a number of features and their associated costs to
every boundary element and nding the minimumcost paths via dynamic
programming When the cursor moves close to the true boundary the live
wire snaps on to the boundary If the result is satisfactory the user deposits
the cursor at which point the location of the cursor becomes the starting
point of a new livewire segment The entire boundary of the ROI is specied
as a set of livewire segments The features used include the image intensity
values on the inside and outside of the boundary several gradient measures
and the distance from the boundary in the preceding slice in the case of seg
mentation of D images Falc%ao et al indicated that the method could assist
in fast and accurate segmentation of ROIs in D and D medical images after
some training
Fusion of multiple results of segmentation The segmentation of re
gions such as masses and tumors in mammograms is complicated by several
factors such as poor contrast superimposition of tissues imaging geometry
and the invasive nature of cancer In such cases it may be desirable to apply
several image processing techniques based upon dierent principles and image
characteristics and then to combine or fuse the multiple results obtained into
a single region or contour It could be expected that the result so obtained
would be better than any of the individual results A fuzzy fusion operator
that can fuse a contour and a fuzzy region to produce another fuzzy region is
described in Section
Application Detection of the Spinal Canal
In an application to analyze CT images of neuroblastoma see Sec
tion the spinal canal was observed to interfere with the segmentation of
the tumor using the fuzzy connectivity algorithm In order to address
this problem a method was developed to detect the center of the spinal canal
in each CT slice grow the D region containing the spinal canal and remove
the structure The initializing seeds for the regiongrowing procedure were
automatically obtained with the following procedure
The outer region in the CT volume containing materials outside the patient
the skin and peripheral fat was rst segmented and removed The
CT volume was then thresholded at HU to detect the highdensity
bone structures All voxels not within mm from the inner boundary of the
peripheral fat layer were rejected Regions were grown using each remaining
voxel and all of the resulting regions were merged to form the bone volume
The inclusion criteria were in terms of the CT values being within
HU with HU being the standard deviation of bone and spatial
connectivity
The resulting CT volume was cropped to limit the scope of further analysis
as follows The width of the image was divided into three equal parts and
the outer thirds were rejected The height of the image was divided into six
equal parts and the lower fourth and fth parts were included in the cropped
region In the interslice direction the rst ! of the slices were removed
and the subsequent slices were included in the cropped volume Figure
b and Figure b show the eect of cropping on two CT slices
The cropped binarized bone volume was subjected to a D derivative oper
ator to produce the edges of the bone structures The vertebral column is not
continuous but made up of interlocking elements As a result the boneedge
map could be sparse Figures c and c show the boneedge maps
of the binarized bone volume related to the corresponding CT images in parts
b of the same gures
The Hough transform for the detection of circles as described in Sec
tion and Equation was applied to each slice of the boneedge map
The radius in the Hough space was limited to the range mm Because
of the possibility of partial structures and edges in a given image the global
maximum in the Hough space may not relate to the inner circular edge of the
spinal canal as desired In order to obtain the center and radius of the ROI
the CT values of bone marrow HU and HU and the spinal
canal HU and HU were used as constraints If the center
of the circle corresponding to the Houghspace maximum was not within the
specied HU range the circle was rejected and the next maximum in the
Hough space was evaluated This process was continued until a suitable circle
was detected
Figure d shows the besttting circle detected drawn on the original
CT image When the bone structure is clearly delineated the besttting
circle approximates well the spinal canal boundary Figure d shows
four circles related to the maximum in the corresponding Hough space and
the subsequent three peaks the related Houghspace slices are shown in Fig
ure The besttting circle which was not given by the global maximum
in the Hough space was obtained by applying the constraints dened above
The centers of the circles detected as above were used as the seed voxels
in a fuzzy connectivity algorithm to segment the spinal canal The mean and
standard deviation required for this procedure were estimated using a
neighborhood around each seed voxel The spinal canal detected over all of
the CT slices available for the case illustrated in Figures and is
shown in Figure The spinal canal volume was then removed from the
CT volume resulting in improved segmentation of the tumor volume
Application Detection of the Breast Boundary in
Mammograms
Identication of the breast boundary is important in order to demarcate the
breast region on a mammogram The inclusion of this preliminary procedure
in CAD systems can avoid useless processing time and data storage By
identifying the boundary of the breast it is possible to remove any artifact
present outside the breast area such as patient markings often highintensity
regions and noise which can aect the performance of image analysis and
pattern recognition techniques Identication and extraction of the eective
breast region is also important in PACS and telemammography systems
The prole of the breast has been used as additional information in dif
ferent tasks in mammography Bick et al and Byng et al for
example used the skinair boundary information to perform density correc
tion of peripheral breast tissue on digital mammograms which is aected
by the compression procedure applied during imaging Chandrasekhar and
Attikiouzel discussed the importance of the skinair boundary prole
as a constraint in searching for the nipple location which is often used as a
reference point for registering mammograms taken at dierent times of the
same subject Other groups have used the breast boundary to perform reg
istration between left and right mammograms in the process of detection of
asymmetry
Most of the works presented in the literature to identify the boundary of
the breast are based upon histogram analysis
which may be critically dependent upon the thresholdselection process and
the noise present in the image Such techniques as discussed by Bick et
al may not be robust for a screening application Ferrari et al
proposed active contour models especially designed to be locally adaptive
for identication of the breast boundary in mammograms The methods
including several preprocessing steps are described in the following sections
a b
c d
FIGURE
a A CT image slice of a patient with neuroblastoma Pixel width
mm The tumor boundary drawn by a radiologist is shown in black
b Cropped area for further processing c Edge map of the binarized bone
volume d Besttting circle as determined by Houghspace analysis The
circle has a radius of mm Figures courtesy of the Alberta Childrens
Hospital and HJ Deglint
a b
c d
FIGURE
a A CT image slice of a patient with neuroblastoma Pixel
width mm b Cropped area for further processing c Edge map
of the binarized bone volume d The four circles corresponding to the rst
four peaks in the Hough space The radii of the circles range from mm
to mm See also Figure Figures courtesy of the Alberta Childrens
Hospital and HJ Deglint
a b
c d
FIGURE
a Edge map of the binarized bone volume related to the ROI of the CT
image in Figure Houghspace slices related to the detection of circles of
radius b c and d pixels with the pixel width being mm
Figures courtesy of HJ Deglint
FIGURE
D rendition of the spinal canal detected for the case related to the CT slices
shown in Figures and The spinal canal is shown in a relatively dark
shade of gray against the bone volume for reference
Detection using the traditional active deformable con
tour model
In an initial study Ferrari et al used the traditional active deformable
contour model or snakes for the detection of the breast boundary The
method summarized in Figure is composed of six main stages as de
scribed in the following paragraphs
Stage The image contrast is enhanced by using a simple logarithmic
operation A contrastcorrection step using a simple logarithmic operation
as
gx y log fx y
is applied to the original image fx y gx y is the transformed image This
operation for dynamic range compression although applied to the whole im
age signicantly enhances the contrast of the regions near the breast bound
ary in mammograms which are characterized by low density and poor de
nition of details The rationale behind the application of this pro
cedure to the image is to determine an approximate breast contour as close
as possible to the true breast boundary The eect of this procedure can be
seen by comparing the original and the enhanced images in Figures a
and b
Stage A binarization procedure using the LloydMax quantization al
gorithm is applied to the image see Section for details Figure
c shows the binarized version of the image in part b of the same
gure
Stage Spurious details generated by the binarization step are removed
by using a morphological opening operator with a circular structuring
element with a diameter of pixels Figures cd show the result of
the binarization procedure for the mammogram in Figure a before and
after the application of the morphological opening operator respectively
Stage After the binarization procedure an approximate contour C
appr
of the breast is extracted by using the chaincode method see Section
for details The starting point of C
appr
is obtained by following the horizontal
path that starts at the centroid of the image and is directed toward the chest
wall until a background pixel is found This procedure avoids selecting an
initial boundary from artifacts or patient labels that may be present in the
image Four control points see Figure e are automatically determined
and used to limit the breast boundary The points are dened as N the
topleft corner pixel of the boundary loop N the farthest point on the
boundary from N in terms of the Euclidean distance through the breast
N the lowest pixel on the lefthand edge of the boundary loop and N the
farthest point on the skinair boundary loop from N
Stage Pixels along lines of length pixels length cm at a sam
pling resolution of m are identied at each point of the approximate
FIGURE
Flowchart of the procedures for identication of the skinair boundary of the
breast in mammograms Reproduced with permission from RJ Ferrari RM
Rangayyan JEL Desautels RA Borges and AF Fr&ere Identication of
the breast boundary in mammograms using active contour models Medical
and Biological Engineering and Computing
c
IFMBE
Figure a
Figure b
Figure c
Figure d
Figure e
Figure f
Figure g
skinair boundary in the original image in the direction normal to the bound
ary see Figure f The graylevel histogram of the pixels along each
normal line is computed and the skinair intersection is dened as the rst
pixel while traversing along the normal line from inside the breast toward
the outside that has the gray level associated with the maximum value in the
histogram as illustrated in Figure This procedure was designed in order
to provide a close estimate to the true skinair boundary see Figure g
and thereby reduce the chances of the active contour used in the next stage
converging to a wrong contour
Stage The traditional snakes model is applied to dene the true breast
boundary The contour determined in the previous stage is used as the input
to a traditional parametric active contour or snakes model The contour
is moved through the spatial domain of the image in order to minimize the
energy functional
E
Z
n
jv
sj
jv
sj
o
E
ext
fvsg
ds
where and are weighting parameters that control respectively the ten
sion and rigidity of the snake The v
s and v
s values denote the rst and
second derivatives of vs with respect to s where vs indicates the continu
ous representation of the contour and s represents distance along the contour
h
FIGURE
Results of each stage of the method for identication of the breast bound
ary a Image mdb from the MiniMIAS database b Image after
the logarithmic operation cd Binary image before and after applying
the binary morphological opening operator e Control points N to N
automatically determined used to limit the breast boundary f Normal
lines computed from each pixel on the skinair boundary g Boundary after
histogrambased analysis of the normal lines h Final boundary Repro
duced with permission from RJ Ferrari RM Rangayyan JEL Desautels
RA Borges and AF Fr&ere Identication of the breast boundary in mam
mograms using active contour models Medical and Biological Engineering
and Computing
c
IFMBE
a
b
FIGURE
a Prole of a sample normal line used to determine an approximate skinair
boundary The symbol indicates the skinair intersection determined in
Stage of the method b Histogram computed from a Reproduced with
permission from RJ Ferrari RM Rangayyan JEL Desautels RA Borges
and AF Fr&ere Identication of the breast boundary in mammograms using
active contour models Medical and Biological Engineering and Computing
c
IFMBE
normalized to the range The external energy function E
ext
fvsg is
derived from the original image fx y as
E
ext
x y krfx yk
where r is the gradient operator The values and were
experimentally derived by Ferrari et al based upon the approximate bound
ary obtained in the previous stage the quality of the external force derived
from the original image and the nal contours obtained
Results of application to mammograms Sixtysix images from the
MiniMIAS database were used to assess the performance of the method
The results were subjectively analyzed by an expert radiologist According
to the opinion of the radiologist the method accurately detected the breast
boundary in images and reasonably well in images The nal contour
for the mammogram in Figure a is given in part h of the same gure
The method failed in ve images because of distortions and artifacts present
near the breast boundary see Figure Limitations of the method exist
mainly in Stages and Although Stage helps in obtaining a good breast
contour it is timeconsuming and may impose a limitation in practical ap
plications The traditional snakes model used in Stage is not robust the
method is not locally adaptive in the presence of noise and artifacts and has
a short range of edge capture
In order to address these limitations Ferrari et al proposed an improved
method described in the following section by replacing the traditional snakes
algorithm with an adaptive active deformable contour model specically de
signed for the application to detect breast boundaries The algorithm includes
a balloon force in an energy formulation that minimizes the in uence of the
initial contour on the convergence of the algorithm In the energy formulation
the external energy is also designed to be locally adaptive
Adaptive active deformable contour model
The improved method to identify the breast boundary is summarized in the
owchart in Figure Stages of the initial method described in
the preceding section are used to nd an approximate breast boundary The
approximate contour V v
v
v
N
with an ordered collection of N
points v
i
x
i
y
i
i N is obtained from Stage of the previous
method by sampling the approximate contour C
appr
see Figure a and
used as the initial contour in the adaptive active deformable contour model
AADCM see Figure b Only ! of the total number of points
present in C
appr
are used in the sampled contour
In the AADCM which combines several characteristics from other known
active contour models the contour is moved through the spatial
domain of the image in order to minimize the following functional of energy
Figure a
E
total
N
X
i
E
internal
v
i
E
external
v
i
where and are weighting parameters that control the internal and external
energies E
internal
and E
external
respectively at each point v
i
The internal
energy is composed of two terms
E
internal
v
i
a E
continuity
v
i
b E
balloon
v
i
This energy component ensures a stable shape for the contour and constrains
to keep constant the distance between the points in the contour In the work
of Ferrari et al the weighting parameters a and b were initially set to unity
a b because the initial contours present smooth shapes and are close
to the true boundary in most cases For each element mn in a neighborhood
of pixels of v
i
the continuity term e
cmn
v
i
is computed as
e
cmn
v
i
lV
jp
mn
v
i
v
i
v
i
j
where lV
N
P
N
i
jv
i
v
i
j
is a normalization factor that makes the
continuity energy independent of the size location and orientation of V
b
FIGURE
Result of the segmentation algorithm showing wrong convergence of the breast
contour to a region of high gradient value a Breast boundary detected su
perimposed on the original image mdb from the MiniMIAS database
b Details of the breast contour attracted to the image identication marker
corresponding to the boxed region in a Compare with Figure Repro
duced with permission from RJ Ferrari RM Rangayyan JEL Desautels
RA Borges and AF Fr&ere Identication of the breast boundary in mam
mograms using active contour models Medical and Biological Engineering
and Computing
c
IFMBE
Figure a
p
mn
v
i
is the point in the image at the position m n in the neigh
borhood of v
i
and cos
N
is a constant factor to keep the location
of the minimum energy lying on the circle connecting v
i
and v
i
in the
case of closed contours see Figure
The balloon force is used to force the expansion of the initial contour toward
the breast boundary In the work of Ferrari et al the balloon force was
made adaptive to the magnitude of the image gradient causing the contour
to expand faster in homogeneous regions and slower near the breast boundary
The balloon energy term e
bmn
v
i
is dened as
e
bmn
v
i
n
i
fv
i
p
mn
v
i
g
where n
i
is the outward unit normal vector of V at the point v
i
and the
symbol indicates the dot product n
i
is computed by rotating the vector
t
i
v
i
v
i
jv
i
v
i
j
v
i
v
i
jv
i
v
i
j
which is the tangent vector at the point v
i
by
o
see Figure
The external energy is based upon the magnitude and direction of the image
gradient and is intended to attract the contour to the breast boundary It is
dened as
e
emn
v
i
n
i
rffp
mn
v
i
g
b
FIGURE
a Approximate breast contour obtained from Stage of the method de
scribed in Section for the image mdb b Sampled breast contour
used as the input to the AADCM Reproduced with permission from RJ
Ferrari RM Rangayyan JEL Desautels RA Borges and AF Fr&ere
Identication of the breast boundary in mammograms using active contour
models Medical and Biological Engineering and Computing
c
IFMBE
FIGURE
Characteristics of the continuity energy component in the adaptive active
deformable contour model Figure adapted with permission from BT Mack
iewich
FIGURE
Characteristics of the balloon energy component in the adaptive active de
formable contour model Figure adapted with permission from BT Mack
iewich
where rffp
mn
v
i
g is the image gradient vector at m n in the
neighborhood of v
i
see Figure The direction of the image gradient is
used to avoid the attraction of the contour by edges that may be located near
the true breast boundary such as identication marks and small artifacts see
Figure In this situation the gradient direction at the position m n
on an edge near the breast boundary and the direction of the unit normal of
the contour will have opposite signs which makes the functional of energy
present a large value at m n
Minimization of the energy functionals In order to allow comparison
between the various energy components described above in the work of Ferrari
et al each energy parameter was scaled to the range as follows
E
continuity
v
i
e
cmn
v
i
e
cmin
v
i
e
cmax
v
i
e
cmin
v
i
E
balloon
v
i
e
bmn
v
i
e
bmin
v
i
e
bmax
v
i
e
bmin
v
i
krfv
i
k
krfk
max
E
external
v
i
e
emn
v
i
e
emin
v
i
maxe
emax
v
i
e
emin
v
i
krfk
max
Here e
min
and e
max
indicate the minimum and maximum of the corresponding
energy component in the neighborhood of v
i
krfk
max
is the maximum
gradient magnitude in the entire image
Ferrari et al used the Greedy algorithm proposed by Williams and Shah
to perform minimization of the functional of energy in Equation
FIGURE
Characteristics of the external energy component in the adaptive active de
formable contour model Figure adapted with permission from BT Mack
iewich
Figure a
Although this algorithm has the drawback of not guaranteeing a global
minimum solution it is faster than the other methods proposed in the litera
ture such as dynamic programming variational calculus and nite elements
It also allows the insertion of hard constraints such as curvature evaluation
as described below
Convergence of the AADCM is achieved in two stages by smoothing the
original image with two dierent Gaussian kernels dened with
x
y
pixels and
x
y
pixels At each stage the iterative process is
stopped when the total energy of the contour increases between consecutive
iterations This coarsetone representation is expected to provide more sta
bility to the contour
In order to allow the deformable contour to adjust to corner regions such
as the upperright limit of the breast boundary a constraint was inserted at
the end of each iteration by Ferrari et al to relax the continuity term dened
in Equation The curvature value Cv
i
at each point v
i
of the contour
was computed as
Cv
i
sin
u
i
ku
i
k
u
i
ku
i
k
where u
i
v
i
v
i
is the vector joining two neighboring contour elements
and is the external angle between such vectors sharing a common contour
b
FIGURE
Application of the gradient direction information to avoid the attraction of the
boundary to objects near the true boundary a Breast boundary detected
automatically superimposed on the original image mdb from the Mini
MIAS database b Details of the detected breast boundary close to the
image identication marker corresponding to the boxed region in the original
image Compare with Figure Reproduced with permission from RJ
Ferrari RM Rangayyan JEL Desautels RA Borges and AF Fr&ere
Identication of the breast boundary in mammograms using active contour
models Medical and Biological Engineering and Computing
c
IFMBE
element This denition of curvature as discussed byWilliams and Shah
has three important advantages over other curvature measures it requires
only simple computation gives coherent values and depends solely on relative
direction
At each v
i
the weight values for the continuity term and the external energy
are set respectively to zero a and to twice the initial value
if Cv
i
Cv
i
and Cv
i
Cv
i
and Cv
i
T The thresh
old value T was set equal to which corresponds to an external angle of
approximately
o
According to Williams and Shah this value of the
threshold has been proven experimentally to be suciently large to dieren
tiate between corners and curved lines Figures b and c illustrate an
example without and with the curvature constraint to correct corner eects
Figure a
In the work of Ferrari et al the weighting parameters and in Equa
tion were initialized to and respectively for each contour element
This set of weights was determined experimentally by using a set of images
randomly selected from the MiniMIAS database not including any im
age in the test set used to evaluate the results A larger weight was given to
the gradient energy to favor contour deformation toward the breast boundary
rather than smoothing due to the internal force
b c
FIGURE
Example of the constraint used in the active contour model to prevent smooth
ing eects at corners a Original image the box indicates the region of
concern b c Details of the breast contour without and with the con
straint for corner correction respectively Reproduced with permission from
RJ Ferrari RM Rangayyan JEL Desautels RA Borges and AF Fr&ere
Identication of the breast boundary in mammograms using active contour
models Medical and Biological Engineering and Computing
c
IFMBE
Results of application to mammograms
Ferrari et al applied their methods to images randomly chosen from the
MiniMIAS database All images were MLO views with m sampling
interval and bit graylevel quantization For reduction of processing time all
images were downsampled with a xed sampling distance so that the original
images corresponding to a matrix size of pixels were transformed
to pixels The results obtained with the downsampled images were
mapped to the original mammograms for subsequent analysis and display
The results were evaluated in consultation with two radiologists experienced
in mammography
The test images were displayed on a computer monitor with a diagonal size
of cm and dot pitch of mm By using the Gimp program the
contrast and brightness of each image were manually enhanced so that the
breast contour could be easily visualized The breast boundary was manually
drawn under the supervision of a radiologist and the results printed on paper
by using a laser printer with dpi resolution The zoom option of the Gimp
program was used to aid in drawing the contours The breast boundaries of
all images were visually checked by a radiologist using the printed images
hardcopy along with the displayed images softcopy the assessment was
recorded for analysis
The segmentation results related to the breast contours detected by image
processing were evaluated based upon the number of falsepositive FP and
falsenegative FN pixels identied and normalized with reference to the cor
responding areas demarcated by the manually drawn contours The reference
area for the breast boundary was dened as the area of the breast image
delimited by the handdrawn breast boundary The FP and FN average per
centages and the corresponding standard deviation values obtained for the
images were ! and ! respectively Thirtythree images
presented both FP and FN percentages less than ! images presented
FP and FN percentages between ! and ! the FP and FN percentages
were greater than ! for images The most common cause of FN pixels
was related to missing the nipple region as illustrated by the example in Fig
ure c By removing Stage used to approximate the initial contour to
the true breast boundary see Figure the average time for processing
an image was reduced from min to min However the number of im
ages where the nipple was not identied increased The AADCM performed
successfully in cases where a small bending deformation of the contour was
required to detect the nipple see Figure
The method for the detection of the breast boundary was used as a pre
processing step in the analysis of bilateral asymmetry by Ferrari et al
see Section The method may also be used in other applications such as
image compression by using only the eective area of the breast and image
registration
Figure a
Figure b
c
FIGURE
Results obtained for the image mdb from the MiniMIAS database
a Original image b Handdrawn boundary superimposed on the histo
gramequalized image c Breast boundary detected superimposed on the
original image Reproduced with permission from RJ Ferrari RM Ran
gayyan JEL Desautels RA Borges and AF Fr&ere Identication of the
breast boundary in mammograms using active contour models Medical and
Biological Engineering and Computing
c
IFMBE
Figure a
Figure b
c
FIGURE
Results obtained for the image mdb from the MiniMIAS database
a Original image b Handdrawn boundary superimposed on the
histogramequalized image c Breast boundary detected superimposed on
the original image Reproduced with permission from RJ Ferrari RM Ran
gayyan JEL Desautels RA Borges and AF Fr&ere Identication of the
breast boundary in mammograms using active contour models Medical and
Biological Engineering and Computing
c
IFMBE
Application Detection of the Pectoral Muscle in
Mammograms
The pectoral muscle represents a predominant density region in most MLO
views of mammograms and can aect the results of image processing meth
ods Intensitybased methods for example can present poor performance
when applied to dierentiate dense structures such as the broglandular disc
or small suspicious masses because the pectoral muscle appears at approxi
mately the same density as the dense tissues of interest in the image The
inclusion of the pectoral muscle in the image data being processed could also
bias the detection procedures Another important need to identify the pec
toral muscle lies in the possibility that the local information of its edge along
with internal analysis of its region could be used to identify the presence
of abnormal axillary lymph nodes which may be the only manifestation of
occult breast carcinoma in some cases
Karssemeijer used the Hough transform and a set of threshold values
applied to the accumulator cells in order to detect the pectoral muscle Ayl
ward et al used a gradientmagnitude ridgetraversal algorithm at small
scale and then resolved the resulting multiple edges via a voting scheme in
order to segment the pectoral muscle Ferrari et al proposed a tech
nique to detect the pectoral muscle based upon the Hough transform
which was a modication of the method proposed by Karssemeijer
However the hypothesis of a straight line for the representation of the pec
toral muscle is not always correct and may impose limitations on subsequent
stages of image analysis Subsequently Ferrari et al proposed another
method based upon directional ltering using Gabor wavelets this method
overcomes the limitation of the straightline representation mentioned above
The Houghtransformbased and Gaborwaveletbased methods are described
in the following sections
Detection using the Hough transform
The initial method to identify the pectoral muscle proposed by Ferrari et
al summarized in the owchart in Figure starts by automatically
identifying an appropriate ROI containing the pectoral muscle as shown in
Figure To begin with an approximate breast contour delimiting the
control points is obtained by using a method for the detection of the skinair
boundary described in Section The six control points NN
used to dene the ROI are determined as follows in order
N the topleft corner pixel of the boundary loop
N the lowest pixel on the lefthand edge of the boundary loop
N the midpoint between N and N
N the farthest point on the boundary from N in terms of the Eu
clidean distance through the breast if this point is not located on the
upper edge of the mammogram it is projected vertically to the upper
edge
N the point that completes a rectangle with N N and N not
necessarily on the boundary loop
N the farthest point on the boundary loop from N In the case of
the mammogram in Figure the points N and N have coincided
The ROI is dened by the rectangular region delimited by the points N
N N and N as illustrated in Figure Although in some cases this
region may not include the total length of the pectoral muscle the portion of
the muscle present is adequate to dene a straight line to represent its edge
By limiting the size of the ROI as described above the bias that could be
introduced by other linear structures that may be present in the broglandular
disc is minimized Diering from the method of Karssemeijer Ferrari et
al did not use any threshold value in order to reduce the number of
unlikely pectoral lines Instead geometric and anatomical constraints were
incorporated into the method as follows
The pectoral muscle is considered to be a straight line limited to an
angle between
o
and
o
with the angle computed as indicated in
Figure Mammograms of right breasts are ipped mirrored before
performing pectoral muscle detection
The pectoral line intercepts the line segment N N as indicated in
Figure
The pectoral line is present in partial or total length in the ROI dened
as the rectangular region delimited by the points N N N and N
as illustrated in Figure
The pectoral muscle appears on mammograms as a dense region with
homogeneous graylevel values
After selecting the ROI a Gaussian lter with
x
y
pixels is used
to smooth the ROI in order to remove the highfrequency noise in the image
The Hough transform is then applied to the Sobel gradient of the ROI to
detect the edge of the pectoral muscle The representation of a straight line
for the Hough transform computation is specied as
x x
cos y y
sin
where x
y
is the origin of the coordinate system dened as the center of
the image and and represent respectively the distance and angle between
FIGURE
Procedure for the identication of the pectoral muscle by using the Hough
transform Reproduced with permission from RJ Ferrari RM Rangayyan
JEL Desautels RA Borges and AF Fr&ere Automatic identication of
the pectoral muscle in mammograms IEEE Transactions on Medical Imag
ing
c
IEEE
Figure a
x
y
and the coordinates x y of the pixel being analyzed as illustrated
in Figure The Hough accumulator is quantized to bins of
o
each
by using the constraint j
xy
j
o
where
xy
is the orientation of the
Sobel gradient at the pixel x y In the work of Ferrari et al
the accumulator cells were incremented using the magnitude of the gradient
instead of unit increments thus pixels with a strong gradient have larger
weights Only values of in the range
o
o
were considered in the
analysis because the pectoral muscle of the mammogram was positioned on
the lefthand side of the image before computing the Hough transform see
Figure
After computing the accumulator cell values in the Hough space a lter
ing procedure is applied to eliminate all lines pairs of parameters and
that are unlikely to represent the pectoral muscle In this procedure all lines
intercepting the top of the image outside the N N line segment see Fig
ure or with slopes outside the range
o
o
are removed In the
present notation the x axis corresponds to
o
and the chest wall is positioned
on the lefthand side see Figure Each remaining accumulator cell is
also multiplied by the factor
A
b
FIGURE
a Image mdb from the MiniMIAS database b Approximate boundary
of the breast along with the automatically determined control points N N
used to limit the ROI rectangle marked for the detection of the pectoral mus
cle Reproduced with permission from RJ Ferrari RM Rangayyan JEL
Desautels RA Borges and AF Fr&ere Automatic identication of the pec
toral muscle in mammograms IEEE Transactions on Medical Imaging
c
IEEE
FIGURE
Coordinate system used to compute the Hough transform The pectoral mus
cle line detected is also shown Reproduced with permission from RJ Ferrari
RM Rangayyan JEL Desautels RA Borges and AF Fr&ere Automatic
identication of the pectoral muscle in mammograms IEEE Transactions
on Medical Imaging
c
IEEE
where and
are respectively the mean and the variance of the graylevel
values in the area A of the pectoral muscle see Figure dened by the
straight line specied by the parameters and This procedure was applied
in order to enhance the Hough transform peaks that dene regions with the
property stated in Item of the list provided earlier in this section page
The weight related to the area was designed to dierentiate the true pectoral
muscle from the pectoralis minor the latter could present a higher contrast
than the former in some cases albeit enclosing a smaller area than the former
Finally the parameters and of the accumulator cell with the maximum
value are taken to represent the pectoral muscle line Figure shows the
Hough accumulator cells at the dierent stages of the procedure described
above for the mammogram in Figure a The pectoral muscle line de
tected for the mammogram in Figure a is shown in Figure
Detection using Gabor wavelets
In an improved method to detect the pectoral muscle edge summarized by the
owchart in Figure Ferrari et al designed a bank of Gabor lters
to enhance the directional piecewiselinear structures that are present in an
ROI containing the pectoral muscle Figure a illustrates a mammogram
from the MiniMIAS database Figure b shows the ROI dened
automatically as the rectangle formed by the chest wall as the lefthand edge
and a vertical line through the uppermost point on the skinair boundary
drawn along the entire height of the mammogram as the righthand edge
to be used for the detection of the pectoral muscle Diering from the
Houghtransformbased method described in Section the ROI here is
dened to contain the entire pectoral muscle region because the straightline
representation is not used Decomposition of the ROI into components with
dierent scale and orientation is performed by convolution of the ROI image
with a bank of tunable Gabor lters The magnitude and phase components of
the ltered images are then combined and used as input to a postprocessing
stage as described in the following paragraphs
Gabor wavelets A D Gabor function is a Gaussian modulated by a
complex sinusoid which can be specied by the frequency of the sinusoid W
and the standard deviations
x
and
y
of the Gaussian envelope as
x y
x
y
exp
x
x
y
y
j Wx
See also Sections and Gabor wavelets are obtained by dilation
and rotation of x y as in Equation by using the generating function
a b c
FIGURE
Hough accumulator cells obtained at three stages of the procedure to detect
the pectoral muscle The contrast of the images has been modied for im
proved visualization a Accumulator cells obtained by using the constraint
j
xy
j
o
and
o
o
b After removing the lines intercepting
the top of the image outside the region dened by the control points N N
see Figure c After applying the multiplicative factor
A
Reproduced with permission from RJ Ferrari RM Rangayyan JEL De
sautels RA Borges and AF Fr&ere Automatic identication of the pectoral
muscle in mammograms IEEE Transactions on Medical Imaging
c
IEEE
FIGURE
Flowchart of the procedure for the identication of the pectoral muscle by
using Gabor wavelets Reproduced with permission from RJ Ferrari RM
Rangayyan JEL Desautels RA Borges and AF Fr&ere Automatic iden
tication of the pectoral muscle in mammograms IEEE Transactions on
Medical Imaging
c
IEEE
a b
FIGURE
a Image mdb from the MiniMIAS database b The ROI used to search
for the pectoral muscle region dened by the chest wall and the upper limit
of the skinair boundary The box drawn in a is not related to the ROI in
b Reproduced with permission from RJ Ferrari RM Rangayyan JEL
Desautels RA Borges and AF Fr&ere Automatic identication of the pec
toral muscle in mammograms IEEE Transactions on Medical Imaging
c
IEEE
mn
x y a
m
x
y
a m n integers
x
a
m
x x
cos
n
y y
sin
n
y
a
m
x x
sin
n
y y
cos
n
where x
y
is the center of the lter in the spatial domain
n
n
K
n K K is the number of orientations desired andm and n indicate
the scale and orientation respectively The Gabor lter used by Ferrari et
al is expressed in the frequency domain as
u v
u
v
exp
uW
u
v
v
where
u
x
and
v
y
The design strategy used by Ferrari et
al was to project the lters so as to ensure that the halfpeak magnitude
supports of the lter responses in the frequency spectrum touch one another
as shown in Figure By doing this it can be ensured that the lters
will capture the maximum information with minimum redundancy In order
for the designed bank of Gabor lters to be a family of admissible D Ga
bor wavelets the lters x y must satisfy the admissibility condition
of nite energy which implies that their Fourier transforms are pure
bandpass functions having zero response at DC This condition is achieved by
setting the DC gain of each lter as which causes the lters not
to respond to regions with constant intensity
The following formulas provide the lter parameters
u
and
v
a
U
h
U
l
S
u
a U
h
a
p
ln
v
tan
K
h
U
h
u
U
h
ln
i
h
ln
ln
u
U
h
i
where U
l
and U
h
denote the lower and upper center frequencies of interest
The K and S parameters are respectively the number of orientations and
the number of scales in the desired multiresolution decomposition procedure
The sinusoid frequency W is set equal to U
h
and m S
In the application being considered interest lies only in image analysis
without the requirement of exact reconstruction or synthesis of the image from
FIGURE
Bank of Gabor lters designed in the frequency domain Each ellipse repre
sents the range of the corresponding lter response from to in squared
magnitude only one half of the response is shown for each lter The sam
pling of the frequency spectrum can be adjusted by changing the U
l
U
h
S
and K parameters of the Gabor wavelets Only the lters shown shaded are
used to enhance the directional piecewiselinear structures present in the ROI
images The frequency axes are normalized Reproduced with permission
from RJ Ferrari RM Rangayyan JEL Desautels RA Borges and AF
Fr&ere Automatic identication of the pectoral muscle in mammograms
IEEE Transactions on Medical Imaging
c
IEEE
the ltered components Therefore instead of using the wavelet coecients
Ferrari et al used the magnitude of the lter response computed as
a
mn
x y
fx y
even
mn
x y
where
even
mn
x y indicates the evensymmetric part of the complex Gabor
lter fx y is the ROI being ltered and represents D convolution The
phase and magnitude images indicating the local orientation were composed
by vector summation of the K ltered images chapter see also
Figure
The area of each ellipse indicated in Figure represents the frequency
spectrum covered by the corresponding Gabor lter Once the range of the
frequency spectrum is adjusted the choice of the number of scales and ori
entations is made in order to cover the range of the spectrum as required
The choice of the number of scales S and orientations K used by Ferrari
et al for detecting the pectoral muscle was based upon the resolution re
quired for detecting oriented information with high selectivity The
spatialfrequency bandwidths of the simple and complex cells in mammalian
visual systems have been found to range from to octaves clustering
around octaves and octaves and their angular bandwidth is expected
to be smaller than
o
By selecting U
l
U
h
S
and K for processing mammographic images Ferrari et al
indirectly adjusted the Gabor wavelets to have a frequency bandwidth of ap
proximately one octave and angular bandwidth of
o
In the work of Ferrari et al all images were initially oriented so that the
chest wall was always positioned on the lefthand side then the pectoral
muscle edge in correctly acquired MLO views will be located between
o
and
o
p Here the orientation of the pectoral muscle edge is
dened as the angle between the horizontal line and an imaginary straight
line representing the pectoral muscle edge For this reason Ferrari et al
used only the Gabor lters with the mean orientation of their responses in
the image domain at
o
o
and
o
the corresponding frequencydomain
responses are shown shaded in Figure
Postprocessing and pectoral muscle edge detection In the method
of Ferrari et al after computing the phase and magnitude images by vector
summation the relevant edges in the ROI are detected by using an algorithm
proposed by Ma and Manjunath for edge ow propagation as described
below The magnitude Ax y and phase x y at each image location x y
are used to represent the edge ow vector instead of using a predictive coding
model as initially proposed by Ma and Manjunath The phase at each point
in the image is propagated until it reaches a location where two opposite
directions of ow encounter each other as follows
Set n and E
x y Ax y cosx y Ax y sinx y
Set the edge ow vector E
n
x y at iteration n to zero
At each image location x y identify the neighbor x
y
that has the
same direction as that of the edge ow vector E
n
x y The direction
is computed as tan
y
y
x
x
If E
n
x
y
E
n
x y
then E
n
x
y
E
n
x
y
E
n
x y
else E
n
x y E
n
x y E
n
x y
where the symbol indicates the dotproduct operation
If nothing has been changed
then stop iterating
else go to Step and repeat the procedure
Figures b and c illustrate an example of the orientation map before
and after applying the edge ow propagation procedure to the image shown
in part a of the same gure
a b c
FIGURE
a Region indicated by the box in Figure a containing a part of the
pectoral muscle edge b and c Edge ow map before and after propaga
tion Each arrowhead represents the direction of the edge ow vector at the
corresponding position in the image Reproduced with permission from RJ
Ferrari RM Rangayyan JEL Desautels RA Borges and AF Fr&ere
Automatic identication of the pectoral muscle in mammograms IEEE
Transactions on Medical Imaging
c
IEEE
After propagating the edge ow vector the boundary candidates for the
pectoral muscle are obtained by identifying the locations that have nonzero
edge ow arriving from two opposite directions Weak edges are eliminated by
thresholding the ROI image with a threshold value of ! of the maximum
graylevel value in the ROI Figure shows images resulting from each
stage of the procedure for pectoral muscle detection
In order to connect the disjoint boundary segments that are usually present
in the image after the edge ow propagation step a halfelliptical neighbor
hood is dened with its center located at each boundary pixel being processed
The halfellipse is adjusted to be proportional to the length of the contour
with R
C
len
and R
pixels where R
and R
are respectively
the major and minor axes of the halfellipse and C
len
is the length of the
boundary with its major axis oriented along the direction of the contour
line If an ending or a starting pixel of an unconnected line segment is found
in the dened neighborhood it is connected by linear interpolation The iter
ative method stops when all disjoint lines are connected this procedure took
to iterations with the images used in the work of Ferrari et al Next
the false edges that may result in the ltered images due to structures inside
the broglandular disc or due to the ltering process see Figure cd
are removed by checking either if their limiting points are far away from the
upper and lefthand side of the ROI or if the straight line having the same
slope as the pectoral muscle edge candidate intercepts outside the upper and
lefthand limits of the ROI Finally the largest line in the ROI is selected to
represent the pectoral muscle edge see Figure e
Results of application to mammograms
Ferrari et al tested their methods with images randomly chosen from the
MiniMIAS database All images were MLO views with m sam
pling interval and bit graylevel quantization For reduction of processing
time all images were downsampled with a xed sampling distance so that the
original images with the matrix size of pixels were transformed
to pixels The results obtained with the downsampled images were
mapped to the original mammograms for subsequent analysis
and display
The results obtained with both of the Houghtransformbased and Gabor
waveletbased methods were evaluated in consultation with two radiologists
experienced in mammography The test images were displayed on a computer
monitor with diagonal size of cm and dot pitch of mm By using the
Gimp program the contrast and brightness of each image were manually
enhanced so that the pectoral muscle edge could be easily visualized Then
the pectoral muscle edges were manually drawn and the results printed on
paper by using a laser printer with dpi resolution The zoom option of
the Gimp program was used to aid in drawing the edges The pectoral muscle
edges of all images were checked by a radiologist using the printed images
a b
Figure c d
e
FIGURE
Result of each stage of the Gaborwaveletbased method a ROI used
b Image magnitude after ltering and vector summation enhanced by
gamma correction cd Results before and after the post
processing stage e Final boundary Reproduced with permission from RJ
Ferrari RM Rangayyan JEL Desautels RA Borges and AF Fr&ere
Automatic identication of the pectoral muscle in mammograms IEEE
Transactions on Medical Imaging
c
IEEE
hardcopy along with the displayed images softcopy the assessment was
recorded for further analysis
The segmentation results were evaluated based upon the number of FP and
FN pixels normalized with reference to the corresponding numbers of pixels
in the regions demarcated by the manually drawn edges The reference region
for the pectoral muscle was dened as the region contained between the left
hand edge of the image and the handdrawn pectoral muscle edge An FP
pixel was dened as a pixel outside the reference region that was included
in the pectoral region segmented An FN pixel was dened as a pixel in the
reference region that was not present within the segmented region Table
provides a summary of the results
TABLE
Average Falsepositive and Falsenegative Rates in the
Detection of the Pectoral Muscle by the houghtransformbased
and Gaborwaveletbased Methods
Method Hough Gabor
FP ! !