Detection of Regions of Interest

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 ! !