Analysis of Oriented Patterns

Analysis of Oriented Patterns

Many images are composed of piecewise linear objects Linear or oriented

objects possess directional coherence that can be quantied and examined to

assess the underlying pattern An area that is closely related to directional im

age processing is texture identication and segmentation For example given

an image of a human face a method for texture segmentation would attempt

to separate the region consisting of hair from the region with skin as well as

other regions such as the eyes that have a texture that is dierent from that

of either the skin or hair In texture segmentation a common approach for

identifying the diering regions is via nding the dominant orientation of the

dierent texture elements and then segmenting the image using this informa

tion The subject matter of this chapter is more focused and concerned with

issues of whether there is coherent structure in regions such as the hair or skin

To put it simply the question is whether the hair is combed or not and if it is

not the degree of disorder is of interest which we shall attempt to quantify

Directional analysis is useful in the eective identication segmentation and

characterization of oriented or weakly ordered texture

Oriented Patterns in Images

In most cases of natural materials strength is derived from highly coherent

oriented bers an example of such structure is found in ligaments

Normal healthy ligaments are composed of bundles of collagen brils that

are coherently oriented along the long axis of the ligament see Figure a

Injured and healing ligaments on the other hand contain scabs of scar ma

terial that are not aligned Thus the determination of the relative disorder

of collagen brils could provide a direct indicator of the health strength and

functional integrity or lack thereof of a ligament similar

patterns exist in other biological tissues such as bones muscle bers and

blood vessels in ligaments as well

Examples of oriented patterns in biomedical images include the following

Fibers in muscles and ligaments see Figure

Fibroglandular tissue ligaments and ducts in the breast see Figures

and

Vascular networks in ligaments lungs and the heart see Figures

and

Bronchial trees in the lungs see Figure

Several more examples are presented in the sections to follow

In manmade materials such as paper and textiles strength usually relies

upon the individual bers uniformly knotting together Thus the strength

of the material is directly related to the organization of the individual bril

strands

Oriented patterns have been found to bear signicant information in several

other applications of imaging and image processing In geophysics the accu

rate interpretation of seismic soundings or stacks is dependent upon the

elimination of selected linear segments from the stacks primarily the ground

roll or lowfrequency component of a seismic sounding Tho

rarinsson et al used directional analysis to discover linear anomalies in

magnetic maps that represent tectonic features

In robotics and computer vision the detection of the objects in the vicinity

and the determination of their orientation relative to the robot are important

in order for the machine to function in a nonstandard environment

By using visual cues in images such as the dominant orientation of a

scene robots may be enabled to identify basic directions such as up and down

Information related to orientation has been used in remote sensing to ana

lyze satellite maps for the detection of anomalies in map data

Underlying structures of the earth are commonly identied by

directional patterns in satellite images for example ancient river beds

Identifying directional patterns in remotely sensed images helps geologists to

understand the underlying processes in the earth that are in action

Because manmade structures also tend to have strong linear segments di

rectional features can help in the identication of buildings roads and urban

features

Images commonly have sharp edges that make them nonstationary Edges

render image coding and compression techniques such as LP coding and

DPCM see Chapter less ecient By dividing the frequency space into di

rectional bands that contain the directional image components in each band

and then coding the bands separately higher rates of compression may be

obtained In this manner directional l

tering can be useful in other applications of image processing such as data

compression

Measures of Directional Distribution

Mardia pointed out that the statistical measures that are commonly used

for the analysis of data points in rectangular coordinate systems may lead to

improper results if applied to circular or directional data Because we do not

usually consider directional components in images to be directed elements or

vectors there should be no need to dierentiate between components that

are at angles and

o

therefore we could limit our analysis to the

semicircular space of

o

o

or

o

o

The rose diagram

The rose diagram is a graphical representation of directional data Corre

sponding to each angular interval or bin a sector a petal of the rose is

plotted with its apex at the origin In common practice the radius of the

sector is made proportional to the area of the image components directed in

the corresponding angle band

The area of each sector in a rose diagram as above varies in proportion

to the square of the directional data In order to make the areas of the

sectors directly proportional to the orientation data the square roots of the

data elements could be related to the radii of the sectors Linear histograms

conserve areas and are comparatively simple to construct however they lack

the strong visual association with directionality that is obtained through the

use of rose diagrams Several examples of rose diagrams are provided in the

sections to follow

The principal axis

The spatial moments of an image may be used to determine its principal axis

which could be helpful in nding the dominant angle of directional alignment

The moment of inertia of an image fx y is at its minimum when the moment

is taken about the centroid x y of the image The moment of inertia of the

image about the line y y cos x x sin passing through x y and

having the slope tan is given by

m

Z

x

Z

y

x x sin y y cos

fx y dx dy

In order to make m

independent of the choice of the coordinates the

centroid of the image could be used as the origin Then x and y

and Equation becomes

m

Z

x

Z

y

x sin y cos

fx y dx dy

m

sin

m

sin cos m

cos

where m

pq

is the p q

th

moment of the image given by

m

pq

Z

x

Z

y

x

p

y

q

fx y dx dy

By denition the moment of inertia about the principal axis is at its minimum

Dierentiating Equation with respect to and equating the result to zero

gives

m

sin m

cos m

sin

or

tan

m

m

m

By solving this equation we can nd the slope or the direction of the principal

axis of the given image

If the input image consists of directional components along an angle only

then If there are a number of directional components at dierent

angles then represents their weighted average direction Evidently this

method cannot detect the existence of components in various angle bands

and is thus inapplicable for the analysis of multiple directional components

Also this method cannot quantify the directional components in various angle

bands

Angular moments

The angular moment M

k

of order k of an angular distribution is dened as

M

k

N

X

n

k

n pn

where n represents the center of the n

th

angle band in degrees pn repre

sents the normalized weight or probability of the data in the n

th

band and N

is the number of angle bands If we are interested in determining the disper

sion of the angular data about their principal axis the moments may be taken

with respect to the centroidal angle M

of the distribution Because the

secondorder moment is at its minimum when taken about the centroid we

could choose k for statistical analysis of angular distributions Hence the

second central moment M

may be dened as

M

N

X

n

n

pn

The use of M

as a measure of angular dispersion has a drawback because

the moment is calculated using the product of the square of the angular dis

tance and the weight of the distribution even a small component at a large

angular distance from the centroidal angle could result in a high value forM

See also Section

Distance measures

The directional distribution obtained by a particular method for an image

may be represented by a vector p

p

p

p

N

T

where p

n

represents the distribution in the n

th

angle band The true distribution of the

image if known may be represented by another vector p

Then the Eu

clidean distance between the distribution obtained by the directional analysis

method p

and the true distribution of the image p

is given as

kp

p

k

v

u

u

t

N

X

n

p

n p

n

This distance measure may be used to compare the accuracies of dierent

methods of directional analysis

Another distance measure that is commonly used is the Manhattan dis

tance dened as

jp

p

j

N

X

n

jp

n p

nj

The distance measures dened above may also be used to compare the

directional distribution of one image with that of another

Entropy

The concept of entropy from information theory see Section can be

eectively applied to directional data If we take pn as the directional PDF

of an image in the n

th

angle band the entropy H of the distribution is given

by

H

N

X

n

pn log

pn

Entropy provides a useful measure of the scatter of the directional elements

in an image If the image is composed of directional elements with a uniform

distribution maximal scatter the entropy is at its maximum if however

the image is composed of directional elements oriented at a single angle or in

a narrow angle band the entropy is close to zero Thus entropy while not

giving the angle band of primary orientation or the principal axis could give

a good indication of the directional spread or scatter of an image

See Figure

Other approaches that have been followed by researchers for the character

ization of directional distributions are numerical and statistical characteri

zation of directional strength morphological operations using a rotating

structural element laser smallangle light scattering and

optical diraction and Fourier analysis

Directional Filtering

Methods based upon the Fourier transform have dominated the area of direc

tional image processing The Fourier transform of an

oriented linear segment is a sinc function oriented in the direction orthogonal

to that of the original segment in the spatial domain see Figure Based

upon this property we can design lters to select linear components at specic

angles However a diculty in using the Fourier domain for directional lter

ing lies in the development of highquality lters that are able to select linear

components without the undesirable eects of ringing in the spatial domain

Schiller et al showed that the human eye contains orientationselective

structures This motivated research on human vision by Marr who

showed that the orientation of linear segments primarily edges is impor

tant in forming the primal sketch Several researchers including Kass and

Witkin Zucker and Low and Coggins used oriented band

pass lters in an eort to simulate the human visual systems ability to identify

oriented structures in images Allen et al developed a verylargescale

integrated VLSI circuit implementation of an orientationspecic retina

Several researchers have used many types of simple lters

with wide passbands at various angles to obtain a redundant decomposition

or representation of the given image Such representations were used to de

rive directional properties of the image For example Kass and Witkin

formed a map of ow lines in the given image and under conformal mapping

obtained a transformation to regularize the ow lines onto a grid The re

sulting transformation was used as a parameter representing the texture of

the image In this manner various types of texture could be recognized or

generated by using the conformal map specic to the texture

Chaudhuri et al used a set of bandpass lters to obtain directional com

ponents in SEM images of ligaments however the lter used was relatively

simple see Sections and Generating highly selective lters in D

is not trivial and considerable research has been directed toward nding gen

eral rules for the formation of D lters Bigun et al developed rules for

the generation of leastsquares optimal beam lters in multiple dimensions

Bruton et al developed a method for designing highquality fan lters

using methods from circuit theory This method results in D recursive lters

that have high directional selectivity and good rollo characteristics and is

described in Section

a b

c d

FIGURE

a A test image with a linear feature b Logmagnitude Fourier spectrum

of the test image in a c Another test image with a linear feature at a

dierent angle d Logmagnitude Fourier spectrum of the test image in b

See also Figure

Sector ltering in the Fourier domain

Fourierdomain techniques are popular methods for directional quantication

of images

The results of research on biological visual systems provide

a biological base for directional analysis of images using lterbased methods

The Fourier transform is the most straightforward method for identifying

linear components The Fourier transform of a line segment is a sinc function

oriented at radians with respect to the direction of the line segment in

the spatial domain see Figure This fact allows the selective ltering of

line segments at a specic orientation by ltering the transformed image with

a bandpass lter

Consider a line segment of orientation slope a and yaxis intercept b in

the x y plane with the spatial limits XX and Y Y In order to

obtain the Fourier transform of the image we could evaluate a line integral

in D along the line y ax b To simplify the procedure let us assume

that the integration occurs over a square region with X Y Because the

function fx y is a constant along the line the term fx y in the Fourier

integral can be normalized to unity giving the equation fx y along the

line y ax b Making the substitution x y ba we have the Fourier

transform of the line image given by

F u v

jaj

Z

Y

Y

Z

Y

Y

exp

j

u

y b

a

v y

dy dy

Y

jaj

exp

j

bu

a

sinc

h

u

a

v

Y

i

From the result above we can see that for the image of a line the Fourier

transform is a sinc function with an argument that is a linear combination

of the two frequency variables u v and with a slope that is the negative

reciprocal of the slope of the original line The intercept is translated into

a phase shift of ba in the u variable Thus the Fourier transform of the

line is a sinc function oriented at

o

to the original line centered about the

origin in the frequency domain regardless of the intercept of the original line

This allows us to form lters to select lines solely on the basis of orientation

and regardless of the location in the space domain Spatial components in

a certain angle band may thus be obtained by applying a bandpass lter in

an angle band perpendicular to the band of interest and applying the inverse

transform If we include a spatial oset in the above calculation it would

only result in a phase shift the magnitude spectrum would remain the same

Figure illustrates the ideal form of the fan lter that may used to select

oriented segments in the Fourier domain

FIGURE

Ideal fan lter in the Fourier domain to select linear components oriented

between

o

and

o

in the image plane Black represents the stopband

and white represents the passband The origin u v is at the center

of the gure

Prior to the availability of highspeed digital processing systems attempts

at directional ltering used optical processing in the Fourier domain Ar

senault et al used optical bandpass lters to selectively lter con

tour lines in aeromagnetic maps Using optical technology Duvernoy and

ChalasinskaMacukow developed a directional sampling method to an

alyze images the method involved integrating along an angle band of the

Fouriertransformed image to obtain the directional content This method

was used by DziedzicGoclawska et al to identify directional content in

bone tissue images The need for specialized equipment and precise instru

mentation limits the applicability of optical processing The essential idea of

ltering selected angle bands however remains valid as a processing tool and

is the basis of Fourierdomain techniques

The main problem with Fourierdomain techniques is that the lters do

not behave well with occluded components or at junctions of linear compo

nents smearing of the line segments occurs leading to inaccurate results when

inverse transformed to the space domain Another problem lies in the trun

cation artifacts and spectral leakage that can exist when ltering digitally

which leads to ringing in the inversetransformed image Ringing artifacts

may be avoided by eective lter design but this in turn could limit the

spatial angle band to be ltered

Considerable research has been reported in the eld of multidimensional

signal processing to optimize directionselective lters

Bigun et al addressed the problem of detection of ori

entation in a leastsquares sense with FIR bandpass lters This method has

the added benet of being easily implementable in the space domain Bruton

et al proposed guidelines for the design of stable IIR fan lters Hou

and Vogel developed a novel method of using the DCT for directional

ltering this method uses the fact that the DCT divides the spectrum into

an upper band and a lower band In D such bandsplitting divides the

frequency plane into directional lter bands By selecting coecients of the

DCT the desired spectral components can be obtained Because the DCT

has excellent spectral reconstruction qualities this results in highquality di

rectionally selective lters A limitation of this technique is that the method

only detects the bounding edges of the directional components because the

bandsplitting in the DCT domain does not include the DC component of the

directional elements

In the method developed by Chaudhuri et al a simple decomposition

of the spectral domain into equal angle bands was employed at

o

per

angle band Each sector lter in this design is a combination of an ideal

fan lter a Butterworth bandpass lter a rampshaped lowpass lter and a

raised cosine window as follows

Hf

r

f

r

f

L

f

r

p

f

r

f

H

q

cos

o

B

where

slope of the weighting function

f

r

normalized radial frequency

p

u

v

p order of the highpass lter

q order of the lowpass lter

f

H

upper cuto frequency normalized

f

L

lower cuto frequency normalized

angle of the Fourier transform sample atanvu

o

central angle of the desired angle band

B angular bandwidth and

weighting factor

The combined lter with

o

and B

o

is illustrated in Figure

Filtering an image with sector lters as above results in component

images Each component image contains the linear components of the original

image in the corresponding angle band

Although the directional lter was designed to minimize spectral leakage

some ringing artifacts were observed in the results To minimize the artifacts

a thresholding method was applied to accentuate the linear features in the

image Otsus thresholding algorithm see Section was applied in

the study of collagen ber images by Chaudhuri et al

FIGURE

Directional sector lter in the Fourier domain The brightness is propor

tional to the gain Figure courtesy of WA Rolston

Thresholding of the component images

Many methods are available for thresholding images with an optimal thresh

old for a given application The component images that result

from the sector lters as described in Section possess histograms that

are smeared mainly due to the strong DC component that is present in most

images Even with highquality lters the DC component appears as a con

stant in all of the component images due to its isotropic nature This could

pose problems in obtaining an eective threshold to select linear image fea

tures from the component images On the other hand the removal of the DC

component would lead to the detection of edges and the loss of information

related to the thickness of the oriented patterns

Otsus method of threshold selection is based upon discriminant mea

sures derived from the graylevel PDF of the given image Discriminant crite

ria are designed so as to maximize the separation of two classes of pixels into

a foreground the desired objects and a background

Consider the graylevel PDF pl of an image with L gray levels l

L If the PDF is divided into two classes C

and C

sepa

rated by a threshold k then the probability of occurrence

i

of the class C

i

i f g is given by

k P C

k

X

l

pl k

k P C

L

X

lk

pl k

and the class mean levels

i

for C

i

i f g are given by

k

k

X

l

l P ljC

k

X

l

l

pl

k

k

k

and

k

L

X

lk

l P ljC

L

X

lk

l

pl

k

T

k

k

where

k

k

X

l

pl

and

k

k

X

l

l pl

are the cumulative probability and rstorder moment of the PDF pl up to

the threshold level k and

T

L

X

l

l pl

is the average gray level of the image

The class variances are given by

k

k

X

l

l

k

P ljC

k

X

l

l

k

pl

k

and

k

L

X

lk

l

k

P ljC

L

X

lk

l

k

pl

k

Using the discriminant criterion

B

k

T

where

B

k

k

k

T

k

k

T

and

T

L

X

l

l

T

pl

Otsus algorithm aims to nd the threshold level k that maximizes the dis

criminant criterion given in Equation Maximizing reduces to max

imizing

B

because the value

T

does not vary with the threshold value k

The optimal threshold value k

is given as

k

arg

max

kL

B

k

Otsus method of thresholding performs well in binarizing a large class of

images

Example Chaudhuri et al applied the directional ltering procedure

described above to a test pattern with line segments of various length width

and gray level at four dierent angles namely

o

o

o

and

o

as

shown in Figure a Signicant overlap was included in order to test

the performance of the ltering procedures under nonideal conditions The

logmagnitude Fourier spectrum of the test image is shown in part b of the

gure directional concentrations of energy are evident in the spectrum The

component image obtained using the ltering procedure for the angle band

o

o

in the Fourier domain is shown in Figure c It is evident

that only those lines oriented at

o

in the image plane have been passed

along with some artifacts The corresponding binarized image using the

threshold value given by Otsus method described above is shown in part d

of the gure Parts e and f of the gure show the binarized components

extracted from the test image for the angle bands

o

o

and

o

o

in the image plane

A close inspection of the component images in Figure indicates that

regions of overlap of lines oriented at dierent directions contribute to each

direction The

o

component image in Figure f has the largest error

due to its low gray level and the large extent of overlap with the other lines

The areas of the line segments extracted by the ltering procedure had errors

with respect to the known areas in the original test image of

and for the

o

o

o

and

o

components respectively

The results of application of methods as above for the directional analysis of

collagen bers in ligaments are described in Section

Design of fan lters

Fan lters are peculiar to D ltering there is no direct D analog of this

type of ltering This fact presents some diculties in designing fan lters

because we cannot easily extend the wellestablished concepts that are used to

design D lters The main problem in the design of fan lters lies in forming

the lter at the origin u v or the DC point in the Fourier domain

At the DC point the ideal fan lter structure has a knife edge which makes

the lter nonanalytic that is if one were to approach the origin in the Fourier

domain from any point within the passband of the fan the limit would ideally

a b

c d

Figure e f

FIGURE

a A test image with overlapping directional components at

o

o

o

and

o

b Logmagnitude Fourier spectrum of the test image Results of

directional ltering with the angle bands specied in the image domain

c

o

o

d Result in c after thresholding and binarization e

o

o

binarized f

o

o

binarized Reproduced with permission from

S Chaudhuri H Nguyen RM Rangayyan S Walsh and CB Frank A

Fourier domain directional ltering method for analysis of collagen alignment

in ligaments IEEE Transactions on Biomedical Engineering

c

IEEE

be unity On the other hand the limit as one approaches the origin from a

point within the stopband should be zero

Various methods have been applied to overcome the problem stated above

such as the use of spectrumshaping smoothing functions with the ideal fan

lter for example Chaudhuri et al used the Butterworth and raised

cosine functions as described in Section However the performance

of even the best spectrumshaping function is limited by the tight spectral

constraints imposed by the nonanalytic point at the origin To obtain better

spectral shaping as in D highorder FIR lters or loworder IIR lters may

be used

The discontinuity at the origin u v is the main problem in the

design of recursive fan lters With nonrecursive or FIR lters this problem

does not result in instability With IIR lters instability can occur if the lters

are not properly designed Stability of lters is usually dened as bounded

input boundedoutput BIBO stability Filters that are BIBO stable

ensure that all inputs that are not innite in magnitude will result in outputs

that are bounded

D lters are commonly derived from real rational continuous functions

of the form

T s

s

Qs

s

P s

s

P

M

m

P

N

n

q

mn

s

m

s

n

P

M

m

P

N

n

p

mn

s

m

s

n

where s

and s

are the Laplace variables the function T s

s

is the Laplace

transformed version of the D partial dierential equation that is related to

the required lter response Qs

s

is the numerator polynomial resulting

from the Laplace transform of the forward dierential forms expressed as a

sum of products in s

and s

with the coecients q

mn

M

and N

represent

the order of the polynomial Q in m and n respectively P s

s

is the de

nominator polynomial obtained from the Laplace transform of the backward

dierential forms expressed as a sum of products in s

and s

with the coef

cients p

mn

and M

and N

represent the order of the polynomial P in m

and n respectively The corresponding frequency response function T u v is

obtained by the substitution of s

j u and s

j v

The discontinuous requirement in the continuous prototype lter at the

origin results in the lter transfer function T s

s

having a nonessential

singularity of the second kind at the origin A nonessential singularity of the

second kind occurs when the numerator and the denominator polynomials

P s

s

and Qs

s

in Equation approach zero at the same frequency

location a

a

resulting in T a

a

The discrete form of the function in Equation is obtained through the

D version of the bilinear transform in D given as

s

i

z

i

z

i

for i

to obtain the following discrete version of the lter

Hz

z

Bz

z

Az

z

P

M

m

P

N

n

b

mn

z

m

z

n

P

M

m

P

N

n

a

mn

z

m

z

n

where the orders of the polynomials M

N

M

and N

are dierent from

the corresponding limits of the continuousdomain lter in Equation due

to the bilinear transform

Filter design using nonessential singularities The D lter design

method of Bruton and Bartley views the nonessential singularity in

herent to fan lters not as an obstacle in the design process but as being

necessary in order to generate useful magnitude responses The method relies

on classical electrical circuit theory and views the input image as a surface

of electrical potential The surface of electrical potential is then acted upon

by a D network of electrical components such as capacitors inductors and

resistors the components act as integrators dierentiators and dissipators

respectively The central idea is to construct a network of components that

will not add energy to the input that is to make a completely passive cir

cuit The passiveness of the resulting circuit will result in no energy being

added to the system that is the lter is nonenergic and thus will ensure

that the lter is stable Bruton and Bartley showed that the necessary

condition for a lter to be stable is that the admittance matrix that links

the current and voltage surfaces must have negative Toeplitz symmetry with

reactive elements supplied by inductive elements that satisfy the nonenergic

constraint

The nonenergic condition ensures that the lter is stable because it implies

that the lter is not adding any energy to the system The maximum amount

of energy that is output from the lter is the maximum amount put into the

system by the input image The derivation given by Bruton and Bartley

is the starting point of a numerical method for designing stable recursive

D lters The derivation shows that recursive lters can be built reliably as

long as the condition above on the admittance matrix is met

Bruton and Bartley provided the design and coecients of a narrow

o

fanstop lter obtained using a numerical optimization method with the

TABLE

Coecients of the Discretedomain Fan Filter with a

o

Fan Stopband

b

mn

n n n

m

m

m

m

m

m

a

mn

n n n

m

m

m

m

m

m

Data courtesy of NR Bartley

condition described above added The method results in a recursive lter

of small order with remarkable characteristics A lter of fth order in z

and second order in z

was designed using this method The corresponding

coecients of the discrete function Hz

z

as in Equation are listed

in Table The coecients in the numerator and denominator each add

up to zero at z

and z

conrming that the lter conforms to the

requirement of the knifeedge discontinuity

Rotation of the lter and image The fan lter design algorithm of

Bruton and Bartley provides lters only for a specic angle band in

the above case for a

o

bandstop lter centered at

o

in the Fourier domain

In order to obtain lters with dierent central orientations it is necessary to

perform a rotation of the prototype lter This may be achieved by using the

following substitution in the analog prototype lter transfer function

s

s

cos s

sin

s

s

where is the amount of rotation desired The discrete version of the lter

is obtained through the bilinear transform The rotation step above is not

the usual rotational transformation for lters but it is necessary to use this

transformation in order to ensure that the lter is stable If the normal

rotational transformation were to be used s

would also be rotated as

s

s

sin s

cos

Then values of s

could turn out to be negative this would indicate that there

would be energy added to the system which would make the lter unstable

Suppose that the prototype lter of the form shown in Equation given

by T

s

s

and with the corresponding frequency response function given

by T

u v is bounded by the straight lines L

and L

passing through the

origin at angles of

p

and

p

with the central line of the lter CL

o

where u as shown in Figure a The lines L

and L

are given by

u cos

p

v sin

p

L

u cos

p

v sin

p

L

As a result of the transformation in Equation the center of the passband

of the rotated frequency response T

r

u v is given as T

u

v

T

u cos

c

v sin

c

v Similarly the straight lines L

and L

are rotated to the straight

lines given by

u cos

p

cos

c

v sin

c

cos

p

sin

p

L

u cos

p

cos

c

v sin

c

cos

p

sin

p

L

see Figure b

A limitation to lter rotation as above is that rotating the lter by more

that

o

would result in a loss of symmetry about the central line of the lter

The rotational warping eect may be compensated for in the prototype lter

T

s

s

In the work of Rolston the prototype lter was rotated by

o

in either direction to obtain lters covering an angle band of

o

o

o

and

o

o

in the Fourier domain Filtering in the range

o

o

was achieved by rotating the image by

o

before passing it through the same

lters as above

The fan lter as above has unit gain which has some drawbacks as well

as some advantages An advantage is that features in the ltered component

images have no more intensity than the corresponding features in the original

image so that image components that exist in the original image in a particu

lar angle band will not be attenuated at all or will be attenuated only slightly

This also limits the number of iterations necessary to attenuate outofband

components for a certain class of images The unit gain is an advantage with

images that have only a small depth of eld because a global threshold can be

used for all component images to obtain a representation of the components

that exist in each specic angle band

a

b

FIGURE

a Original fan lter b The fan lter after rotation by the transformation

given in Equation Figure courtesy of WA Rolston

Example A test image including rectangular patches of varying thickness

and length at

o

o

o

and

o

with signicant overlap is shown in

Figure a The fan lter as described above for the central angle of

o

was applied to the test image Because the fan lter is a fanstop lter

the output of the lter was subtracted from the original image to obtain

the fanpass response As evident in Figure b the other directional

components are still present in the result after one pass of the lter due to

the nonideal attenuation characteristics of the lter The lter however can

provide improved results with multiple passes or iterations The result of

ltering the test image nine times is shown in Figure c The result

possesses good contrast between the desired objects and the background and

may be thresholded to reject the remaining artifacts

Gabor Filters

Most of the directional fan and sector lters that have been used in the

Fourier domain to extract directional elements are not analytic functions

This implies that lter design methods in D are not applicable in D Such

lters tend to possess poor spectral response and yield images with not only

the desired directional elements but also artifacts

One of the fundamental problems with Fourier methods of directional lter

ing is the diculty in resolving directional content at the DC point the origin

in the Fourier domain The design of highquality fan lters requires

a

b c

FIGURE

a A test image with overlapping directional components at

o

o

o

and

o

Results of fan ltering at

o

after b one pass c nine passes Figure

courtesy of WA Rolston

conicting constraints at the DC point approaching the DC point from any

location within the passband requires the lter gain to converge to unity how

ever approaching the same point from any location in the stopband requires

the lter gain to approach zero In terms of complex analysis this means

that the lter is not analytic or that it does not satisfy the CauchyRiemann

equations This prevents extending results in D to problems in D

The Gabor lter provides a solution to the problem mentioned above by

increasing the resolution at the DC point Gabor lters are complex sinu

soidally modulated Gaussian functions that have optimal localization in both

the frequency and space domains Gabor lters have been used for tex

ture segmentation and discrimination and

may yield better results than simple Fourier methods for directional ltering

Timelimited or spacelimited functions have Fourier spectra of unlimited

extent For example the timelimited rectangular pulse function transforms

into the innitely long sinc function On the other hand the timeunlimited

sine function transforms into a delta function with innite resolution in the

Fourier domain Innitely long functions cannot be represented in nite cal

culating machinery Gabor suggested the use of timelimited functions

as the kernels of a transform instead of the unlimited sine and cosine functions

that are the kernel functions of the Fourier transform The functional nature

of the Fourier transform implies that there exists an uncertainty principle

similar but not identical to the wellknown Heisenberg uncertainty principle

of quantum mechanics Gabor showed that complex sinusoidally modulated

Gaussian basis functions satisfy the lower bound on the fundamental uncer

tainty principle that governs the resolution in time and frequency given by

t f

where t and f are time and frequency resolution respectively The un

certainty principle implies that there is a resolution limit between the spatial

and the Fourier domains Gabor proved that there are functions that can

form the kernel of a transform that exactly satisfy the uncertainty relation

ship The functions named after Gabor are Gaussianwindowed sine and co

sine functions By limiting the kernel functions of the Fourier transform with

a Gaussian windowing function it becomes possible to achieve the optimal

resolution limit in both the frequency and time domains The size of the

Gaussian window function needs to be used as a new parameter in addition

to the frequency of the sine and cosine functions

Gabor functions provide optimal joint resolution in both the Fourier and

time domains in D and form a complete basis set through phase shift and

scaling or dilation of the original mother basis function The set of functions

forms a multiresolution basis that is commonly referred to as a wavelet basis

formalized by Mallat

Daugman extended Gabor functions to D as D sinusoidal plane

waves of some frequency and orientation within a D Gaussian envelope Ga

bor functions have also been found to provide good models for the receptive

elds of simple cells in the striate cortex for this reason there has

been a signicant amount of research conducted on using the functions for

texture segmentation analysis and discrimination

The extension of the principle above to D leads to spacelimited plane

waves or complex exponentials Such an analysis was performed by Daug

man The uncertainty relationship in D is given by

x y u v

where x and y represent the spatial resolution and u and v represent

the frequency resolution The D Gabor functions are given as

hx y gx

y

expj Ux V y

x

y

x cos y sinx sin y cos

where x

y

are the x y coordinates rotated by an arbitrary angle

gx y

exp

x

y

is a Gaussian shaping window with the aspect ratio and U and V are the

center frequencies in the u v frequency plane An example of the real part

of a Gabor kernel function is given in Figure with U

V and with reference to Equations and Another

Gabor kernel function is shown in gray scale in Figure

The imaginary component of the Gabor function is the Hilbert transform

of its real component The Hilbert transform shifts the phase of the original

function by

o

resulting in an odd version of the function

The Gabor transform is not a transform as such that is there is usually

no transform domain into which the image is transformed The frequency

domain is usually divided into a symmetric set of slightly overlapping regions

at octave intervals Examples of the ranges related to a few Gabor functions

are shown in Figure see also Figures and It is evident

that Gabor functions act as bandpass lters with directional selectivity

Multiresolution signal decomposition

Multiresolution signal analysis is performed using a single prototype function

called a wavelet Fine temporal or spatial analysis is performed with con

tracted versions of the wavelet on the other hand ne frequency analysis is

performed with dilated versions The denition of a wavelet is exible and

requires only that the function have a bandpass transform thus a wavelet at

a particular resolution acts as a bandpass lter The bandpass lters must

FIGURE

An example of the Gabor kernel with U V and

with reference to Equations and Figure courtesy of WA

Rolston

FIGURE

An example of a Gabor kernel displayed as an image Figure courtesy of

WA Rolston

FIGURE

Division of the frequency domain by Gabor lters Two sets of oval regions

are shown in black corresponding to the passbands of three lters in each

set oriented at

o

and

o

In each case the three regions correspond to

three scales of the Gabor wavelets There is a

o

shift between the angles

of corresponding lter functions in the space and frequency domains Figure

courtesy of WA Rolston

have constant relative bandwidth or constant quality factor The importance

of constant relative bandwidth of perceptual processes such as the auditory

and visual systems has long been recognized Multiresolution analysis

has also been used in computer vision for tasks such as segmentation and

object recognition The analysis of nonstationary signals

often involves a compromise between how well transitions or discontinuities

can be located and how nely longterm behavior can be identied This

is reected in the abovementioned uncertainty principle as established by

Gabor

Gabor originally suggested his kernel function to be used over bandlimited

equally spaced areas of the frequency domain or equivalently with constant

window functions This is commonly referred to as the shorttime Fourier

transform STFT for shorttime analysis of nonstationary signals

The D equivalent of the STFT is given by

F

S

x

y

u v

Z

x

Z

y

fx y wx x

y y

expj ux vy dx dy

where w is a windowing function and f is the signal image to be analyzed

The advantage of shorttime or movingwindow analysis is that if the energy

of the signal is localized in a particular part of the signal it is also localized

to a part of the resultant D space x

y

u v The disadvantage of this

method is that the same window is used at all frequencies and hence the

resolution is the same at all locations in the resultant space The uncertainty

principle does not allow for arbitrary resolution in both of the space and

frequency domains thus with this method of analysis if the window function

is small the largescale behavior of the signal is lost whereas if the window

is large rapid discontinuities are washed out In order to identify the ne or

smallscale discontinuities in signals one would need to use basis functions

that are small in spatial extent whereas functions of large spatial extent

would be required to obtain ne frequency analysis By varying the window

function one will be able to identify both the discontinuous and stationary

characteristics of a signal The notion of scale is introduced when the size

of the window is increased by an order of magnitude Such a multiresolution

or multiscale view of signal analysis is the essence of the wavelet transform

Wavelet decomposition in comparison to STFT analysis is performed over

regions in the frequency domain of constant relative bandwidth as opposed to

a constant bandwidth

In the problem of determining the directional nature of an image we have

the discontinuity in the frequency domain at the origin or DC to overcome

Wavelet analysis is usually applied to identify discontinuities in the spatial

domain however there is a duality in wavelet analysis provided by the un

certainty principle that allows discontinuity analysis in the frequency domain

as well In order to analyze the discontinuity at DC largescale or dilated

versions of the wavelet need to be used This is the dual of using contracted

versions of the wavelet to analyze spatial discontinuities

The wavelet basis is given by

h

x

y

x y

p

h

x x

y y

where x

y

and

are real numbers and h is the basic or mother wavelet

For large values of

and

the basis function becomes a stretched or ex

panded version of the prototype wavelet or a lowfrequency function whereas

for small

and

the basis function becomes a contracted wavelet that is

a short highfrequency function

The wavelet transform is then dened as

F

W

x

y

p

Z

x

Z

y

fx y

h

x x

y y

dx dy

From this denition we can see that wavelet analysis of a signal consists of

the contraction dilation and translation of the basic mother wavelet and

computing the projections of the resulting wavelets on to the given signal

Formation of the Gabor lter bank

In the method proposed by Bovik et al the given image is convolved

with the complex Gabor kernel and the maximum magnitude of the result

is taken as an indicator to identify changes in the dominant orientation of

the image In the work of Rolston and Rangayyan this method

was observed to fail in the presence of broad directional components The

real component of the Gabor lter acts as a matched lter to detect broad

directional components and thus is better suited to the identication of such

regions

The parameters of Gabor lters that may be varied are as follows With

reference to Equations and the parameter species the spatial

extent of the lter species the aspect ratio of the lter that modulates the

value If the parameter in Equation need not be specied

because gx y is then isotropic In the frequency domain such a lter results

in an oriented lter occupying the middle subsection of the corresponding ideal

fan lter with the orientation being specied by tan

VU see Figure

These parameters will then completely specify the Gabor lter bank

In the directional analysis algorithm proposed by Rolston and Rangayyan

only the real component of the Gabor wavelet is used with

and the primary orientation given by tan

VU

o

o

o

and

o

A given image is analyzed by convolving bandlimited and

decimated versions of the image with the same analyzing wavelet

When a decimated image is convolved with a lter of constant spatial ex

tent relative to the original image the lter is eectively scaled larger with

respect to the decimated image The advantage of this procedure is that

lters with larger values or with center frequencies closer to DC can be

simulated instead of resorting to using lters of larger spatial extent Filters

with larger values correspond to portions of the frequency domain closer

to the DC point This eect is shown in Figure The frequency plane is

completely covered by the decimation and ltering operation Each black oval

in Figure represents the frequency band or region that is being ltered by

each decimation and ltering operation The largest black oval at each orien

tation corresponds to onetoone ltering and the smaller ovals closer to the

origin correspond to higher orders of decimation and ltering Higher levels

of decimation and ltering geometrically approach the DC point Theoreti

cally arbitrary resolution of the DC point can be achieved using this method

however the size of the original image imposes a limiting factor because an

image that has been digitized for example to pixels can only be

decimated a few times before the details of interest are lost

Another advantage of this method is that because ltering is performed

in the spatial domain stable lters are obtained at DC whereas with strict

Fouriertransformbased methods the resolution at DC is reduced to one fre

quency increment or less thereby making the Fourierdomain lters sensitive

Because the lter bank works on decimated images the computational load

of convolution reduces geometrically at each successive stage of decimation

Reconstruction of the Gabor lter bank output

In the directional ltering and analysis procedures proposed by Rolston and

Rangayyan the given image is decimated and convolved at each of

three scales with a lter of xed size Decimation and ltering at each scale

results in equal energy across all of the scales due to the selection of the lter

coecients Thus after interpolation of the decimated and convolved images

the responses at the dierent scales can be added without scaling to obtain

the overall response of the lter at the dierent scales

After obtaining the responses to the lters at

o

o

o

and

o

a vector

summation of the lter responses is performed as shown in Figure The

vector summation is performed at each pixel in the original image domain to

obtain a magnitude and phase or angle at each pixel

FIGURE

Vector summation of the responses of Gabor lters at

o

o

o

and

o

Figure courtesy of WA Rolston

The Gabor lters as described above do not have a perfect reconstruction

condition This results in a small amount of outofband energy interfering

with the reconstruction translated into artifacts in the reconstructed image

A thresholding operation can eectively remove such artifacts Rolston and

Rangayyan set the threshold as the maximum of the sum of the

mean and standard deviation values across all of the component images

Example A test image including rectangular patches of varying thickness

and length at

o

o

o

and

o

with signicant overlap is shown in Figure

a The output at each stage of the Gabor lter as described above for

o

is shown in same gure It is evident that with onetoone ltering the narrow

horizontal elements have been successfully extracted however only the edges

of the broader components are present in the result As the decimation ratio

is increased the broader components are extracted which indicates that the

ltering is eective in lowfrequency bands closer to the DC point

See Sections and for further discussions and results related

to Gabor lters

Directional Analysis via Multiscale Edge Detection

Methods for edge detection via multiscale analysis using LoG functions are

described in Section Liu et al applied further steps to the

edge stability map obtained by this method see Figure to detect linear

segments corresponding to collagen bers in SEM images of ligaments which

are described in the following paragraphs

Estimating the area of directional segments Directional analysis

requires the estimation of the area covered by linear segments in specied

angle bands For this purpose the pattern boundaries obtained by the relative

stability index see Equation may be used for computing the area of

coverage The directional information of a pattern is given by the directions

of the gradients along the detected pattern boundaries

Figure a depicts the approach of Liu et al to area compu

tation where two patterncovered regions are denoted by R

A

and R

B

The

arrows along the boundaries indicate the directions of the gradients which are

computed from the original image on a discretized grid depending upon the

application The use of gradients enables the denition of the region enclosed

by the boundaries It is seen from Figure a that a linear segment can

be identied by a pair of line segments running in opposite directions In

order to identify the region Liu et al proposed a piecewise labeling

procedure that includes two steps line labeling and region labeling In the

linelabeling procedure the full plane is sectioned into eight sectors see

Figure and a set of templates is dened for pixel classication The

relative stability index is scanned left to right and top to bottom To each

element in the relative stability index a line label is assigned according to its

match with one of the templates A structure array is constructed to store

the descriptions of the lines at both pixel and line levels The structure array

contains several description elds such as the line starting location xs ys

the ending location xe ye the orientation and a corner label which is

also a structure array containing the corner location and the lines that form

the corner

a

b c

d e

FIGURE

a A test image with overlapping directional components at

o

o

o

and

o

Results of Gabor ltering at

o

after decimation at b onetoone

c twotoone and d fourtoone e Overall response at

o

after vec

tor summation as illustrated in Figure Figure courtesy of WA Rol

ston

FIGURE

a Computation of the area covered by directional segments The arrows

perpendicular to the pattern boundaries represent gradient directions used for

detecting the interior of the linear segment over which the area is computed

The directional information associated with the pattern is also stored for anal

ysis b Computation of occluded segments based upon the detected Tjoints

The subscripts denote dierent regions and the superscripts denote the line

numbers Reproduced with permission from ZQ Liu RM Rangayyan and

CB Frank Directional analysis of images in scalespace IEEE Transac

tions on Pattern Analysis and Machine Intelligence

c

IEEE

Once the line segments have been labeled a set of region descriptors is gen

erated which includes paired line labels their starting and ending locations

orientation and the area of the region see Figure a In region label

ing a line for example Line

A

is paired with an adjacent line for example

Line

A

having a direction that is in the sector opposite to that of Line

A

see Figure The area of the linear segment R

A

is then computed by

counting the number of pixels contained by the pair of line segments The

orientation of the linear segment is indicated by the orientation of the pair of

line segments For instance if Line

A

and Line

A

form a pair their associated

region descriptor can be dened as

RfA xs ys xe ye

xs ys xe ye

g

where the subscripts and represent Line

A

and Line

A

respectively and

is the area computed for the region R

A

see Figure a

FIGURE

The image plane is divided into eight sectors Line

and Line

form a

pair Reproduced with permission from ZQ Liu RM Rangayyan and

CB Frank Directional analysis of images in scalespace IEEE Transac

tions on Pattern Analysis and Machine Intelligence

c

IEEE

Detection of occluded linear segments It is often the case in natural

images particularly in SEM images of ligaments that as linear patterns in

tersect some segments of a linear pattern will be occluded Analysis based

upon incomplete patterns will introduce errors It is desirable that the oc

cluded linear segments are detected and included in the nal analysis In the

methods of Liu et al a simple interpolation method is used for this

purpose

Occluded segments typically appear as Tjunctions in an edge image See

Chen et al for a discussion on various types of junctions in images with

oriented patterns As described above a corner structure array is generated

along with the line structure array Tjunctions can be readily detected by

inspecting the corners and if necessary linking lines according to the following

procedure The lines that form Tjunctions with a common line see Figure

b are considered to be occluded line segments and are stored in a

Tjunction array structure

Tfk Line

A

Line

A

Line

B

Line

B

Line

k

g

where k indicates the k

th

Tjunction structure and the subscript indicates

the region associated with the common line After all the Tjunction struc

tures are constructed they are paired by bringing together the Tjunction

structures with Line

k

that share the same region Corresponding line ele

ments in paired Tjunction structures are then compared to detect lines that

cut across the common region This is performed by verifying if a line in

one of the Tjunction structures of the pair lies within a narrow coneshaped

neighborhood of the corresponding line in the other Tjunction structure of

the pair If such a line pair is detected across a pair of Tjunction structures

the lines are considered to be parts of a single line with an occluded part under

the common region Furthermore if two such occluded lines form two regions

on either side of the common region the two regions are merged by adding

the occluded region and relabeled as a single region With reference to the

situation depicted in Figure b the above procedure would merge the

regions labeled as R

D

and R

E

into one region including the area occluded in

between them by R

Overview of the algorithm for directional analysis In summary the

method proposed by Liu et al for directional analysis via multiscale

ltering with LoG functions see Section consists of the following steps

Generate a set of zerocrossing maps images

Classify or authenticate the zero crossings

Generate the adjusted zerocrossing maps from the original zerocrossing

maps

Generate a stability map from the set of adjusted zerocrossing maps

Generate the relative stability index map from the stability map

Compute the edge orientation from the relative stability index map and

the original image

Compute the orientational distribution of the segments identied

Compute statistical measures to quantify the angular distribution of

the linear patterns such as entropy and moments as described in Sec

tion

The methods described above were tested with the image in Figure a

The areas of the line segments extracted by the procedures had errors with

respect to the known areas in the original test image of

and for the

o

o

o

and

o

components respectively

Liu et al applied the procedures described above for the analysis

of collagen remodeling in ligaments this application is described in detail in

Section

HoughRadon Transform Analysis

The Hough transform is a method of transforming an image into a parameter

domain where it is easier to obtain the desired information in the image see

Section The main drawback of the Hough transform is that it is primar

ily applicable to binary images as a consequence the results are dependent

upon the binarization method used for segmenting the image Rangayyan and

Rolston proposed the use of a combination of the Hough transform

and the Radon transform see Section that overcomes this drawback

their methods and results obtained are described in the following paragraphs

Limitations of the Hough transform

With reference to Figure we see that a straight line can be specied in

terms of its orientation with respect to the x axis and its distance from

the origin In this form of parameterization any straight line is bounded

in angular orientation by the interval and bounded by the Euclidean

distance to the farthest point of the image from the center of the image The

equation for an arbitrary straightline segment in the image plane is given by

x cos y sin

For a specic point in the image domain x

i

y

i

we obtain a sinusoidal curve

in the Hough domain Each point x

i

y

i

lying on a straight line with

and

in the image domain corresponds to a sinusoidal curve in

the domain specied by

x

i

cos

y

i

sin

Through Equation it is evident that for each point in the image do

main the Hough transform performs a onetomany mapping resulting in a

modulated sum of sinusoids in the Hough domain

The Hough transform is often referred to as a voting procedure where each

point in the image casts votes for all parameter combinations that could have

produced the point All of the sinusoids resulting from the mapping of a

straight line in the image domain have a common point of intersection at

in the Hough domain Linear segments in the spatial domain corre

spond to largevalued points in the Hough domain see Figures and

Thus the problem of determining the directional content of an image becomes

a problem of peak detection in the Hough parameter space

FIGURE

Parameters in the representation of a straight line in the Hough transform

and a ray in the Radon transform Reproduced with permission from RM

Rangayyan and WA Rolston Directional image analysis with the Hough

and Radon transforms Journal of the Indian Institute of Science

c

Indian Institute of Science

From the properties listed above the Hough transform appears to be the

ideal tool for detecting linear components in images However there are some

limitations to this approach The results are sensitive to the quantization

intervals used for the angle and the distance Decreasing the quantization

step for increases the computation time because the calculation for needs

to be performed across each value of and each pixel Another problem with

this method is the crosstalk between multiple straight lines in the Hough

domain If the image contains several lines parallel to the x axis they would

correspond to several peak values in the Hough domain at diering values

for

o

However the Hough transform would also detect false linear

segments for

o

which would show up as smaller peaks at a continuum

of values in the Hough domain see Figure This is caused by the fact

that the Hough transform nds line segments at specic values that are not

necessarily contiguous Another form of crosstalk can occur within a broad

directional element several straight lines may be perceived within a broad

element with angles spread about the dominant orientation of the element as

well as at several other angles see Figure

FIGURE

Crosstalk between multiple lines causing the Hough transform to detect false

lines In the case illustrated several short segments of vertical lines are de

tected in addition to the true horizontal lines

The Hough transform has the desirable feature that it handles the occlusion

of directional components gracefully because the size of the parameter peaks

is directly proportional to the number of matching points of the component

The Hough transform also has the feature that it is robust to the addition

of random pixels from poor segmentation because random image points are

unlikely to contribute coherently to a single point in the parameter space

The Hough and Radon transforms combined

Deans showed that there is a direct relationship between the Hough and

Radon transforms The Hough transform may be viewed as a special case of

the Radon transform but with a dierent transform origin and performed on

a binary image Typically the Radon transform is dened with its transform

origin at the center of the original image whereas the Hough transform is

dened with its transform origin at the location of the image where the row

and column indices are zero Thus the distance as in Equation for a

image for the Hough transform would be calculated relative to the

FIGURE

False detection of straight lines at several angles dashed lines within a broad

linear feature by the Hough transform

point in the original image whereas for the Radon transform the

value would be calculated relative to the point see Figure

In the method proposed by Rangayyan and Rolston a Hough

Radon hybrid transform is computed by updating the

i

i

parameter point

by adding the pixel intensity and not by incrementing by one as with the

Hough transform In this sense brighter lines correspond to larger peaks in

the HoughRadon domain The HoughRadon space is indexed from

o

to

o

along one axis and from N to M

p

N

p

for an image with M rows

and N columns as shown in Figure

The generation of the HoughRadon space produces relative intensities of

the directional features in the given image An example of the HoughRadon

space is shown in Figure for a simple test pattern In directional analysis

it would be of interest to obtain the number of pixels or the percentage of

the image area covered by linear segments within a particular angle band

Therefore it is necessary to form a shadow parameter space with the numbers

of the pixels that are in a particular cell in the parameter space The shadow

parameter space is the Hough transform of the image with no accompanying

threshold

It is necessary to form both the HoughRadon transform space and the

Houghtransform shadow space because performing only the Hough transform

on an unthresholded image will produce for most images a transform with

little information about the image Computing the shadow parameter space

is the same as performing the Hough transform on a thresholded image for

all pixels with a gray level greater than zero The HoughRadon transform

however facilitates the dierentiation between the light and dark regions in an

FIGURE

Mapping of a straight line from the image domain to the HoughRadon space

Reproduced with permission from RM Rangayyan and WA Rolston Di

rectional image analysis with the Hough and Radon transforms Journal of

the Indian Institute of Science

c

Indian Institute of Science

FIGURE

a A test image with ve line segments b The HoughRadon space of

the image c Filtered HoughRadon space d Rose diagram of directional

distribution See also Figure Reproduced with permission from RM

Rangayyan and WA Rolston Directional image analysis with the Hough

and Radon transforms Journal of the Indian Institute of Science

c

Indian Institute of Science

image thus retaining all of the information about the image while performing

the desired transform The Hough transform is needed for further processing

when deriving the numbers of pixels regardless of the related intensity

From the result shown in Figure b we can see the high level of

crosstalk in the upperright quadrant From Figure we see that this sec

tion maps to the angle band

o

o

This is due to the Hough transforms

tendency to identify several lines of varying orientation within a broad linear

segment as illustrated in Figure this is both a strength and a weakness

of the Hough transform A ltering procedure is described in the following

subsection to reduce this eect

Filtering and integrating the Hough Radon space

Although the HoughRadon transform is a powerful method for determining

the directional elements in an image it lacks by itself the means to eliminate

elements that do not contribute coherently to a particular directional pattern

This is due to the transform being performed on all points in the given image

simple integration along a column of the transform space will include all the

points in the image A second step is needed to eliminate those pixels that do

not contribute signicantly to a particular pattern Leavers and Boyce

proposed a simple lter to locate maxima in the Hough space that

correspond to connected collinearities in an edge image space

The lter is derived from the parameterization of lines and the ex

pected shape of the distribution of counts in the accumulator of the Hough

space For a linear element in an image the expected shape is a characteris

tic buttery which is a term commonly used to describe the typical fallo

from a focal accumulator point as shown in Figure It was shown by

Leavers and Boyce that for any line in the image space the extent of

the corresponding buttery in the Hough domain is limited to one radian or

approximately

o

of the corresponding focal accumulator point

The D lter

provides a high positive response to a distribution that has its largest value at

the focal point and falls o to approximately on either side and vanishes

rapidly above and below the focal point A drawback of this lter is that it

was designed for detecting peaks in the Hough space corresponding to lines of

one pixel width In the example shown in Figure b we can see that the

broad directional components in the test image correspond to broad peaks in

the HoughRadon domain This results in the lter of Equation detecting

only the edges of the peaks in the Hough domain an example of this eect is

shown in Figure c

The lter in Equation is also sensitive to quantization of the incre

ments This can be seen in the vertical streaks of intensity in Figure c

The vertical streaks occur at values that correspond to points where the

value of in Equation approaches an integral value Increasing the spa

tial extent of the lter would reduce the sensitivity of the lter to noise as well

as improve the ability of the lter to detect larger components in the Hough

Radon transform space Detecting larger components in the HoughRadon

domain corresponds to detecting broad directional image components

In the method proposed by Rangayyan and Rolston after the

HoughRadon transform has been ltered using the lter in Equation

the result is normalized to the range of to and then multiplied point

bypoint with the shadow Hough transform mentioned earlier This step is

performed in order to obtain the relative strength of the numbers of pixels

at each of the detected peaks This step also reduces the accumulated quan

tization noise from the HoughRadon transformation and the ltering steps

Although peaks may be detected in a region of the HoughRadon domain

there may be few corresponding pixels in the original image that map to such

locations Multiplying noisy peaks by areas that contain few points in the

original image will reduce the nal count of pixels in the corresponding angle

bands

The nal integration step is a simple summation along each of the columns

of the ltered parameter space Because the Hough transform generates a

parameter space that is indexed in the column space from

o

to

o

each of

the columns represents a fraction of a degree depending upon the quantization

interval selected for the transform Also because the Hough transform is a

voting process the peaks selected will contain some percentage of the pixels

that are contained in the directional components

Example Figure shows the results of the methods described above for

a simple test image The problem of the smallextent lter mentioned above

is evident the lter detects only the edges of the transformed components

and neglects the central section of each of these components From the rose

diagram shown in Figure d we can see that the lter has detected

the relative distribution of the linear components however there is a large

amount of smearing of the results leading to poor dierentiation between the

angle bands We can also see from the rose diagram that there is a large

amount of crosstalk around

o

where the result should be zero This eect

is probably due to a combination of crosstalk as well as quantization noise in

the HoughRadon transform

For the ligament image shown in Figure a the method has performed

reasonably well Regardless the results contain artifacts due to the limitation

of the algorithm with broad directional components as described above

See Rangayyan and Krishnan for an application of the HoughRadon

transform for the identication of linear sinusoidal and hyperbolic frequency

modulated components of signals in the timefrequency plane

FIGURE

a An SEM image of a normal ligament with wellaligned collagen bers

b The HoughRadon space of the image c Filtered HoughRadon space

d Rose diagram of directional distribution See also Figure Reproduced

with permission from RM Rangayyan and WA Rolston Directional im

age analysis with the Hough and Radon transforms Journal of the Indian

Institute of Science

c

Indian Institute of Science

Application Analysis of Ligament Healing

The collagenous structure of ligaments Virtually all connective

tissues in the human body have some type of berlled matrix The bers

consist of various sizes and shapes of chemically distinct proteins known as

collagen with augmentation by other brous materials such as elastin

There is a complex interaction between these materials and the nonbrous

ground substance in all tissues water proteoglycans glycoproteins and

glycolipids giving each tissue relatively unique mechanical properties As

with any composite berreinforced material the quantity and quality as

well as the organization of the reinforcing bers have considerable inuence

on the mechanical behavior of ligaments

Ligaments are highly organized connective tissues that stabilize joints Lig

aments normally consist of nearly parallel arrangements of collagen bers that

are attached to bone on both sides of a joint serve to guide the joint through

its normal motions and prevent its surfaces from becoming separated Colla

gen bers and their component brils make up the protenaceous backbone

of ligaments and provide the majority of their resistance to tensile loading

The spatial orientation of collagen brils is an important factor in determining

tissue properties Ligaments need to be loose enough to allow joints to move

but have to be tight enough to prevent the joint surfaces from separating

Injuries to ligaments are common with the normal highly structured tis

sue being replaced by relatively disordered scar tissue The scar tissue in

ligaments has many quantitative and qualitative dierences from the normal

ligament but as with scar in other tissues the relative disorga

nization of its collagenber backbone may be among the most critical The

loose meshwork of the scar may not be able to resist tensile loads within the

same limits of movement and deformation as a normal ligament The injured

or healing joint therefore may become loose or unstable

The ne vascular anatomy of ligaments When a ligament is

damaged by injury the extent of the healing response determines whether and

when normal function of the ligament will return A critical factor thought to

be important for the healing of a ligament is its blood supply which exchanges

oxygen nutrients and proteins with ligament tissue The nature of

ligament vascularity has been qualitatively assessed in a few studies

However despite the potential importance of ligament bloodvessel

anatomy not much quantitative information is available on either normal or

healing vascular anatomy of ligaments

With respect to normal ligament vascularity the medial collateral ligament

MCL of the knee in a rabbit model commonly used for ligament healing

studies has previously been characterized qualitatively Excluding its

bony attachments the MCL complex of the rabbit is composed of two main

tissue types epiligament and ligament tissue proper The epiligament is a

thin layer of loose connective tissue surrounding the supercial surface of the

ligament proper Blood vessels in the normal uninjured ligament tissue

proper appear sparse and are oriented parallel to the long axis of the liga

ment in an organized fashion whereas blood vessels in the normal epiligament

appear more abundant and are oriented in a less organized fashion In

ligament scar tissue early after ligament injury blood vessels have been de

scribed to be larger more abundant and more disorganized early on in the

ligament healing process The need for a greater supply of materials to

the ligament for early healing apparently leads to the formation of many new

blood vessels but with longer term maturation of healing tissue the vascular

supply decreases and vascularity may eventually return to normal

Some of the qualitative and quantitative dierences between normal and

healing ligaments have been described in an animal model by Frank et al

Quantitative studies on collagen ber organization were conducted by Chaud

huri et al Frank et al and Liu et al the methods and results

of these studies are described in Section Eng et al and Bray

et al conducted studies with the aim to develop a method to analyze

quantitatively the variations in vascularity in normal and healing ligaments

to correlate such information with other aspects of ligament healing in a well

characterized ligament healing model to predict ligament healing based upon

the vascular response to injury and to develop better methods to optimize lig

ament vascularity after injury The related methods and results are discussed

in Section

Analysis of collagen remodeling

Tissue preparation and imaging The animal model selected in the studies

of Chaudhuri et al Frank et al and Liu et al was the ruptured

and unrepaired MCL in the sixmonthold female New Zealand white rabbit

Under general anesthesia and with sterile techniques the right MCL was

exposed through longitudinal medial incisions in the skin and fascia The right

MCL was completely ruptured by passing a D braided steel wire beneath the

ligament and failing the ligament with a strong upward pull on both ends of

the suture The left MCL was not ruptured and served as a normal control

The injured MCL was allowed to heal for a specied period The animal

was then sacriced by intravenous injection of mg of phenobarbitol and

the healing right and normal control left MCLs were harvested as follows

The right and left MCLs were exposed through medial incisions in the skin

and fascia The MCLs were xed in situ by dropping a fresh solution of

gluteraldehyde in M cacodylate buer with pH onto their

surfaces The MCLs were then removed at their insertions placed in

gluteraldehyde in M cacodylate buer with pH for three hours

and dehydrated in increasing concentrations of ethanol and

Each xed and dehydrated ligament was then frozen quickly in liquid

nitrogen and fractured longitudinally to expose the internal collagen ber and

component bril arrangement along its length The fractured tissue was then

critically point dried aligned mounted and sprayed with gold! palladium

In each case the longitudinal axis of the tissue was distinguishable at low

magnication so as to allow the orientation of highmagnication photographs

relative to that axis

Specimens were viewed under a Hitachi S SEM In order to select parts

of the ligaments for imaging pairs of x y coordinates were obtained using a

random number generator In every image the vertical axis of the photograph

was aligned with the longitudinal axis of the original ligament tissue A

number of photographs were taken randomly in the midsubstance area of each

of the healing and normal control MCLs at a magnication of This

magnication was experimentally chosen to give a good compromise between

the resolution of the collagen brils and the area of the tissue being sampled

The resulting images were then digitized into arrays

Directional analysis A sample image of a normal ligament is shown in

Figure a Parts b and c of the gure show two binarized component

images obtained via directional ltering for the angle bands

o

o

and

o

o

using the sectorltering methods described in Section Directional

components were obtained over angle bands spanning the full range of

o

o

The fractional bercovered areas in the components are shown in

the form a rose diagram in part d of the gure The rose diagram indicates

that most of the collagen bers in the normal ligament tissue sample are

aligned close to the long axis of the ligament

o

in the image plane

A sample image of a oneweek scar tissue sample is shown in Figure a

It is readily seen that the collagen bers in the scar tissue do not have any

dominant or preferred orientation Two binarized component images for the

angle bands

o

o

and

o

o

are shown in parts b and c of the gure

The rose diagram in part d of the gure shows that the angular distribution

of collagen in the healing tissue is almost uniform or random

Frank et al conducted a detailed assessment of collagen realignment

in healing ligaments in response to three methods of treatment immobiliza

tion of the aected joint for three weeks or six weeks and no immobilization

Scar tissue samples were obtained from three rabbits each at three six and

weeks after injury except in the second case that lacked the threeweek

samples Sample images from the groups with no immobilization and immo

bilization for three weeks are shown in Figure Sets of images were

obtained from each sample at randomly chosen locations Composite rose

diagrams were computed for each group and are shown in Figure for the

groups with no immobilization and immobilization for three weeks Figures

and demonstrate the collagen remodeling or realignment process in

healing ligaments

Plots of the entropy of the rose diagrams for all the cases in the study are

shown in Figure The plots clearly demonstrate a reduction in entropy

indicating a return to orderly structure as the healing time increases Im

mobilization of the aected joint for three weeks after injury has resulted in

a b

c d

FIGURE

a A sample image showing collagen alignment in a normal ligament Bina

rized directional components in the angle band b

o

o

and c

o

o

d Fractional bercovered areas in the form a rose diagram Figure courtesy

of S Chaudhuri

a b

c d

FIGURE

a A sample image showing collagen alignment in ligament scar tissue Bina

rized directional components in the angle band b

o

o

and c

o

o

d Fractional bercovered areas in the form a rose diagram Figure courtesy

of S Chaudhuri

entropy values that are close to the values at weeks in all cases and well

within the range for normal ligaments the shaded region in Figure

The results indicate that immobilization of the aected joint for three weeks

promotes the healing process and that immobilization for the longer period

of six weeks does not provide any further advantage

Among the limitations of the methods described above it should be noted

that there is a certain degree of depth to the micrographs analyzed The con

tribution of bril components at varying depths in the dierent angle bands

is aected by a number of factors the intensity of the electron beam in the

microscope the depth of focus photographic methods and image digitization

parameters The nal thresholding scheme applied to the ltered component

images being adaptive in nature aects the various component images dier

ently depending upon their content this aspect cannot be controlled without

introducing bias Artifacts are also caused by tissue xation and handling

spaces between brils ber disruption and surface irregularities Regard

less the results provide important quantitative information that can assist in

the understanding of ligament structure and healing

Analysis of the microvascular structure

Tissue preparation and imaging In the works of Eng et al and

Bray et al adult New Zealand White rabbits were used Only

monthold females were used to reduce variability between animals due to age

or gender dierences The selected rabbits received a gap injury in the right

MCLs by surgical removal of a mm segment of the tissue see Figure

The remaining gapinjury site was marked by means of four small nylon

sutures attached to the original cut ligament ends The injured animals were

allowed to heal for periods of three six or weeks The MCLs from the

injured animals were then removed as described below Selected uninjured

animals were used as normal controls The right and left MCLs from these

animals were removed by the procedure as described below

The control and healing rabbits were sacriced with an overdose of sodium

pentobarbitol Euthanyl cc kg Using a constant volume pump

Indiaink solution was perfused through the femoral arteries of each hind limb

to vessels in the knee ligaments according to an established protocol

When complete perfusion of the limbs was evident by noting when the nail

beds of claws were blackened perfusion was stopped and the entire hind

limb was removed placed in a container at

o

C for a minimum of four

hours to allow the ink solution to set and subsequently the entire MCL

was removed as shown in Figure The left MCLs of the injured animals

were not considered to be normal due to possible eects of injury to the

opposite contralateral knee and were used as contralateral control MCLs

The ligament scar material that developed over the gapinjury site on the

right MCL was labeled as the midsubstance scar see Figure The two

a b

FIGURE

Sample images showing collagen alignment in ligament samples at three weeks

six weeks and weeks after injury a without immobilization of the af

fected joint b with immobilization of the aected joint for three weeks

Images courtesy of CB Frank See also Figure

a b

FIGURE

Composite rose diagrams showing collagen realignment in ligament samples at

three weeks six weeks and weeks after injury a without immobilization

of the aected joint b with immobilization of the aected joint for three

weeks See also Figure Reproduced with permission from CB Frank

B MacFarlane P Edwards R Rangayyan ZQ Liu S Walsh and R Bray

A quantitative analysis of matrix alignment in ligament scars A comparison

of movement versus immobilization in an immature rabbit model Journal

of Orthopaedic Research

c

Orthopaedic Research

Society

FIGURE

Variation of the entropy of composite rose diagrams with collagen realign

ment in ligament samples at three weeks six weeks and weeks after injury

The vertical bars indicate one standard deviation about the corresponding

means NON without immobilization of the aected joint IMM with

immobilization of the aected joint for three weeks IMM with immobi

lization of the aected joint for six weeks The shaded region indicates the

range of entropy for normal ligament samples See also Figures and

Reproduced with permission from CB Frank B MacFarlane P Edwards

R Rangayyan ZQ Liu S Walsh and R Bray A quantitative analysis of

matrix alignment in ligament scars A comparison of movement versus immo

bilization in an immature rabbit model Journal of Orthopaedic Research

c

Orthopaedic Research Society

ligament sections directly connected to the midsubstance scar were labeled as

the original ligament ends

FIGURE

Gapinjury site in the ligament and the formation of scar A Gap injury cre

ated by removing a mm section of the MCL B Scar after healing C Ex

tracted ligament and its main regions See also Figure Reproduced

with permission from K Eng RM Rangayyan RC Bray CB Frank L

Anscomb and P Veale Quantitative analysis of the ne vascular anatomy

of articular ligaments IEEE Transactions on Biomedical Engineering

c

IEEE

Following removal the ligament samples were xed in paraformalde

hyde frozen in an embedding medium sagitally sectioned at m thick

ness using a Cryostat ReichertJung Frigocut N and contactmounted

maintaining proper orientation This procedure resulted in black inklled

blood vessels in ligament tissue which when viewed under light microscopy

showed black vessels on a clear background of ligament collagen Every sec

tion was divided using a grid system and under a microscope the image of

every tenth eld was photographed Photographs that did not contain visible

blood vessels were discarded The number of discarded blank photographs

was taken into account when estimating the fractional volume of blood vessels

in the ligament The specimen magnication was measured to be and

the scale of the digitized image was m per pixel

Typical images of a normal ligament section and a week healing ligament

section are shown in Figure It is evident that the normal ligament is

relatively avascular and that the blood vessels that exist are aligned along

the length of the ligament the horizontal direction of the image On the

FIGURE

Ligament sectioning procedure for the imaging of vascular anatomy A knee

joint B Extracted ligament and plane of sectioning a MCL complex b Lig

ament c Epiligament d Femur e Tibia f Sectioning imaging plane

See also Figure Reproduced with permission from K Eng RM Ran

gayyan RC Bray CB Frank L Anscomb and P Veale Quantitative

analysis of the ne vascular anatomy of articular ligaments IEEE Transac

tions on Biomedical Engineering

c

IEEE

other hand the scar tissue has a more abundant network of blood vessels to

facilitate the healing process with extensive branching and lack of preferred

orientation

a

b

FIGURE

Microvascular structure in ligaments a normal b week scar Images

courtesy of RC Bray

Binarization of the images Images of ligament sections obtained as

above contained inkperfused blood vessels as well as collagen brils and other

artifacts In order to simplify the image analysis procedure the graylevel im

ages were thresholded to create binary images consisting of only two gray

levels representing the blood vessels and the background The graylevel his

togram for a bloodvessel image was assumed to be bimodal with the rst

peak representing the pixels of blood vessels and the second one represent

ing the background pixels Otsus method see Section for threshold

selection produced binary images with excessive artifacts Threshold selec

tion by histogram concavity analysis a method that locates the locally

signicant minima and maxima in the graylevel histogram of the image and

produces a list of possible thresholds was investigated however it was di

cult to choose the actual threshold to be used from the list Another method

tried was the RutherfordAppleton thresholdselection algorithm which

computes a threshold by using gradient information from the image The best

threshold for the binarization of the bloodvessel images was obtained by using

the RutherfordAppleton algorithm to get a threshold estimate followed by

histogram concavity analysis to ne tune the nal threshold value as follows

The derivatives of the given image fmn were obtained in the x and y

directions as

d

x

mn fmn fmn

and

d

y

mn fm n fm n

The larger of the two derivatives was saved as

dmn maxjd

x

mnj jd

y

mnj

Two sums were computed over the entire image as

S

d

X X

dmn

and

S

df

X X

dmn fmn

The RutherfordAppleton threshold is given as

T

o

S

df

S

d

Another potential threshold was determined by nding the position of maxi

mal histogram concavity A typical graylevel histogram consists of a number

of signicant peaks local maxima and valleys local minima Signicant

peaks may be identied by constructing a convex hull of the histogram which

is dened as the smallest convex polygon

"

hl containing the given histogram

hl where l stands for the graylevel variable The convex hull consists of

straightline segments joining the signicant peaks in the histogram The his

togram concavity at any gray level is dened as the vertical distance between

the convex hull and the histogram that is

"

hlhl Within each straight

line segment of the convex hull the gray level at which the maximal concavity

occurred was labeled as the optimal threshold for that segment

Because the area covered by the blood vessels is small compared to the

area covered by the background in the ligament section images the graylevel

histogram was rst scaled logarithmically to make the histogram peak repre

senting the blood vessels and the background peak closer in height A convex

polygon of the scaled histogram was then constructed The problem of choos

ing between the thresholds of each of the segments of the convex polygon

was addressed by nding a threshold T

o

using the RutherfordAppleton algo

rithm described above The threshold estimate T

o

was found to lie between

the background peak and the peak representing the bloodvessel pixels The

threshold representing the maximal histogram concavity within the convex

hull segment joining these two peaks was chosen to be the threshold value T

c

A threshold was also determined by nding the minimum point in the his

togram between the peaks that represented the bloodvessel and background

pixels This threshold labeled T

m

yielded a smaller value than T

c

because of

the height dierence between the peaks

Comparing the various methods of threshold selection described above it

was observed that T

c

was often too high resulting in an image with artifacts

The threshold T

m

was often too low resulting in the loss of bloodvessel

pixels By using the average of T

c

and T

m

as the nal threshold an acceptable

compromise was reached

A sample image of the vascular structure in a normal ligament is shown

in Figure a The histogram of the image with the various thresholds

described above is shown in Figure The binarized version of the image

is shown in Figure b

Skeletonization Skeletonization makes directional analysis easier by re

ducing the binary bloodvessel patterns to their skeletal patterns with one

pixelthick lines see Section In the work of Eng et al the bina

rized bloodvessel images as in Figure b were reduced to their skeletons

by the method described in Section see Figures and for exam

ples

In order to assist the analysis of both the directionality and the volume of

vascularization an image array containing the diameter of the blood vessel at

each skeleton point was formed and referred to as the diameterproportional

skeleton of the image The diameter at a skeleton point s

i

was obtained as

x y minDs

i

C

where C is the set of contour points of the binary image before skeletonization

and D is the Euclidean distance measure

It is to be noted that during skeletonization a line is shortened by onehalf

of its original thickness at each of its two ends In order to correct for this

during reconstruction and area estimation pixels need to be added to the end

points Because most of the bloodvessel images were observed to have smooth

contours the ends of the line segments were assumed to be semicircular The

areas of such half circles were added at the end points

Directional analysis Skeletonization allows the use of the simple method

of leastsquares linear regression to determine the angle of orientation

a

b

FIGURE

Microvascular structure in a normal ligament sample a original image

b binarized image See Figure for details on the selection of the thresh

old for binarization Reproduced with permission from K Eng RM Ran

gayyan RC Bray CB Frank L Anscomb and P Veale Quantitative

analysis of the ne vascular anatomy of articular ligaments IEEE Transac

tions on Biomedical Engineering

c

IEEE

FIGURE

Logarithmically scaled histogram of the image in Figure a along

with its convex hull and several possible thresholds for binarization RATS

RutherfordAppleton thresholdselection algorithm Reproduced with per

mission from K Eng RM Rangayyan RC Bray CB Frank L Anscomb

and P Veale Quantitative analysis of the ne vascular anatomy of articular

ligaments IEEE Transactions on Biomedical Engineering

c

IEEE

FIGURE

Skeleton of the image in Figure b See also Figure

of each bloodvessel segment in the image In the work of Eng et al

from each point x y in the skeleton image a line segment consisting of

N points was extracted with the center point located at x y If x

i

y

i

i N represent the points in the line segment the slope of the

besttting straight line is given by

m

P

N

i

x

i

P

N

i

y

i

P

N

i

x

i

y

i

h

P

N

i

x

i

i

P

N

i

x

i

It should be noted that when the slope becomes large for a nearly vertical

line segment slope estimation as above becomes inaccurate due to increasing

yaxis errors This error can be obviated by adapting the leastsquares formula

to minimize the xaxis errors if the slope found by Equation is greater

than unity The inverse of the slope is then given by

m

P

N

i

x

i

P

N

i

y

i

P

N

i

x

i

y

i

h

P

N

i

y

i

i

P

N

i

y

i

The angle of the skeleton at the point x y is then given by atanm

The elemental area of the blood vessel at the point x y is

Ax y x yW

where x y is the vessel thickness at x y as given by Equation and

W

cos

if jj

o

sin

if jj

o

The factor W as above in pixels accounts for the fact that diagonally con

nected pixels are farther apart than vertically or horizontally connected pixels

The elemental area was added to the corresponding angle of the histogram

and the process repeated for all points in the skeleton

The overall accuracy of the directional analysis procedure as above was

estimated to be

o

by analyzing various test patterns For this reason

the bloodvessel angular distributions were computed in bins of width

o

Figure shows composite rose diagrams obtained from images from

four normal ligaments and images from three ligament scar samples at

weeks of healing It is evident that the normal ligaments demonstrate a well

contained angular distribution of blood vessels angular SD

o

entropy

out of a maximum of whereas the scar tissues show relatively

widespread distribution angular SD

o

entropy

a b

FIGURE

Angular distributions of blood vessels in a normal ligaments averaged over

images from four ligaments and b week scar tissues from three lig

aments images Reproduced with permission from K Eng RM Ran

gayyan RC Bray CB Frank L Anscomb and P Veale Quantitative

analysis of the ne vascular anatomy of articular ligaments IEEE Transac

tions on Biomedical Engineering

c

IEEE

TABLE

Measures of Entropy and Standard Deviation SD of Composite Angular

Distributions of Blood Vessels in Ligaments

Tissue type Ligaments Images Entropy SD

o

Vasc

NORMAL

Ligament

Epiligament

CONTRALATERAL

Ligament

Epiligament

SCAR

ENDS

Ligament

Epiligament

The maximum possible value for entropy is #SCAR midsubtance scar

#ENDS original ligament ends see Figures and # Vasc per

centage of the analyzed tissue volume covered by the blood vessels detected

Reproduced with permission from K Eng RM Rangayyan RC Bray CB

Frank L Anscomb and P Veale Quantitative analysis of the ne vascular

anatomy of articular ligaments IEEE Transactions on Biomedical Engineer

ing

c

IEEE

In addition to the directional distributions and their statistics entropy and

angular dispersion or standard deviation the relative volume of blood vessels

in the various ligament samples analyzed were computed see Table Using

the twosample Ttest several assertions were arrived at about the relative

volume and organization of blood vessels in normal and healing ligaments see

Table Statistical analysis of the results indicated with condence

that week scars contain a greater volume of blood vessels than normal

ligaments Using entropy as a measure of chaos in the angular distribution of

the bloodvessel segments statistical analysis indicated with condence

that blood vessels in week scars are more chaotic than in normal ligaments

A factor that aects the accuracy in the angular distributions derived as

above is the width of the blood vessels As the thickness of a blood vessel in

creases more material is lost at the ends of the vessels during skeletonization

This loss although corrected for by the addition of semicircular end pieces

could lead to reduced accuracy of the angular distribution Sampling and

quantization errors become signicant when the thickness of blood vessels is

small

TABLE

Results of Statistical Comparison of the Relative Volume of

Vascularization V and the Entropy of the Angular Distribution H of

Various Ligament Samples

Assertion Condence

LIGAMENT

V normal V contralateral

V normal V midsubstance scar

V normal V original ligament ends

V original ligament ends V midsubstance scar

H contralateral H normal

H normal H midsubstance scar

H normal H original ligament ends

H original ligament ends H midsubstance scar

EPILIGAMENT

V normal V contralateral

V normal V original ligament ends

H normal H contralateral

H normal H original ligament ends

Reproduced with permission from K Eng RM Rangayyan RC Bray CB

Frank L Anscomb and P Veale Quantitative analysis of the ne vascular

anatomy of articular ligaments IEEE Transactions on Biomedical Engineer

ing

c

IEEE

It should be observed that blood vessels branch and merge in D within

the ligament The sectioning procedure used to obtain D slices imposes a

limitation in the analysis the segments of the blood vessels that traverse

across the sectioning planes are lost in the procedures for directional analysis

The increased bloodvessel volume and greater directional chaos observed

in scar tissue as compared to normal ligaments indicate that the blood supply

pattern is related to the healing process Interestingly after a week healing

period increases in both the bloodvessel volume and dispersion were also ob

served in the MCL from the opposite knee contralateral as compared to the

uninjured normal MCL The directional chaos of the contralateral

MCL decreased in the ligament region but increased in the epiligament region

of the tissue as compared to the normal tissue This may be attributed to a

possible change in loading conditions on the contralateral limb after injury

or to a nervous response to injury transmitted by neural pathways As

expected both the vascularity and the directional chaos of the blood vessels in

the original ligament ends increased as compared to the normal This shows

that injury to one portion of the tissue has some eect on the vascularity of

connected ligament tissue

Application Detection of Breast Tumors

The detection of breast tumors is dicult due to the nature of mammographic

images and features Many studies have focused on image processing tech

niques for the segmentation of breast parenchyma in order to detect suspicious

masses and reject false positives The detection of masses requires the seg

mentation of all possible suspicious regions which may then be subjected to

a series of tests to eliminate false positives

Early work by Winsberg et al using lowresolution imagery showed

promise in detecting large solitary lesions in mammograms Later on sev

eral studies dealt with the detection of suspicious areas in xeromammograms

Ackerman and Gose developed computer techniques for categorizing sus

picious regions marked by radiologists based upon features in xeromammo

grams Kimme et al proposed an automatic procedure for the detection

of suspicious abnormalities on xeromammograms by identifying breast tissues

and partitioning them into at most sections per image Ten normalized

statistics for each section were used as texture features and the classication

of mammographic sections of eight patients with six bestperforming

features yielded a falsepositive rate of and a falsenegative rate of

Hand et al and Semmlow et al reported on automated meth

ods to detect suspicious regions in xeromammograms The methods included

routines for the detection of the breast boundary and the extraction of suspi

cious areas Classication of normal and abnormal regions was achieved using

global and regional features With a dataset of xeroradiograms from

patients the methods proposed by Hand et al correctly identied

of the suspicious areas presented but resulted in a high falsepositive rate

Lai et al presented a method to detect circumscribed masses They

enhanced the contrast in mammograms using a selective median lter as a

preprocessing step and proposed a method based upon template matching

to identify candidate regions Finally two tests local neighborhood test and

region histogram test were applied to reduce the number of false positives

The method was eective in the detection of specic types of lesions but was

not generally applicable to all types of masses The masses detected included

both benign and malignant types however the results of detection were not

reported separately for the individual mass types

Multiscale methods The features of masses in mammograms vary great

ly in size and shape Several computer techniques have therefore

employed multiscale concepts for the detection of masses in mammographic

images Brzakovic et al proposed fuzzy pyramid linking for mass local

ization and used intensity links of edge pixels in terms of their position at

various levels of resolution of the image relative stability of the links from

one level of resolution to another was assumed They reported detec

tion accuracy with cases containing benign masses malignant tumors and

normal images the numbers of each type were not specied They further

reported to have achieved accuracy in classifying the regions detected as

benign malignant or nontumor tissue by using features based upon tumor

size shape and intensity changes in the extracted regions Details about the

range of the size of the masses detected and their distribution in terms of

circumscribed and spiculated types were not reported

Chen and Lee used multiresolution wavelet analysis and expectation

maximization EM techniques in conjunction with fuzzy Cmeans concepts

to detect tumor edges The method was tested on only ve images and clas

sication of masses as benign or malignant was not performed Li et al

developed a segmentation method based upon a multiresolution Markov ran

dom eld model and used a fuzzy binary decision tree to classify the regions

segmented as normal tissue or mass The algorithm was reported to have

achieved sensitivity with two falsepositive detections per image Us

ing a nonlinear multiscale approach Miller and Ramsey achieved an

accuracy of in detecting malignant tumors in a screening dataset

Qian et al developed three image processing modules for

computerassisted detection of masses in mammograms a treestructured

nonlinear lter for noise suppression a multiorientation directional wavelet

transform and a multiresolution wavelet transform for image enhancement

to improve the segmentation of suspicious areas They performed wavelet

decomposition of the original images and used the lowpasssmoothed images

for detecting the central core portions of spiculated masses Waveletbased

adaptive directional ltering was performed on the highpassdetail images to

enhance the spicules in the mass margins Chang and Laine used coher

ence and orientation measures in a multiscale analysis procedure to enhance

features and provide visual cues to identify lesions

Densitybased methods Masses are typically assumed to be hyperdense

with respect to their surroundings However due to the projection nature of

mammograms overlapped broglandular tissues could also result in high

intensity regions in the image leading to their detection as false positives

Therefore many studies focused on detecting poten

tial densities as an initial step and incorporated modules to reject the false

positives at a later stage Woods and Bowyer detected potential densi

ties by using a local measure of contrast and supplemented the method with a

regiongrowing scheme to nd the extents of the potential densities They fur

ther used a set of features computed from each region in linear and quadratic

neural classiers to classify the regions segmented as masses or false posi

tives In a later work Woods and Bowyer examined ve massdetection

algorithms and described a general framework for detection algorithms com

posed of two steps pixellevel segmentation and regionlevel classication

Their review reveals some fundamental advantages in concentrating eort on

pixellevel analysis for achieving higher detection accuracies

Kok et al and Cerneaz referred to a massdetection approach

that was developed by Guissin and Brady based upon isointensity con

tours and reported that the measures that Guissin and Brady used to reduce

false positives are not suitable for mammography Cerneaz and Brady

proposed methods to remove curvilinear structures including blood ves

sels milk ducts and brous tissues prior to the detection of signicant den

sities in mammograms However the assumption that such a step would not

aect the spicules of tumors is questionable

In the detection scheme proposed by Cerneaz dense regions blobs in

the mammogram were initially segmented by combining the methods proposed

by Guissin and Brady and Lindeberg Novelty analysis methods

were then applied to separate mass regions from the plethora of dense regions

thus segmented Three of the ve features used in the novelty analysis step

included variance in the intensity of a blobs pixels average blob height and

a measure of a blobs saliency A set of images no details were men

tioned about the nature of the distribution of the images from the MIAS

database were used for evaluating the methods after reducing the res

olution to m The detection accuracies were not stated explicitly The

detection approach proposed by Mudigonda et al and described in the

subsections to follow possesses similarities with the method of Guissin and

Brady and Cerneaz and Brady in the initial stage of segmen

tation of isolated regions in the image

Analysis of texture of mass regions Petrick et al reported on

the use of a twostage adaptive densityweighted contrast enhancement lter

in conjunction with a LoG edge detector for the detection of masses In their

approach the original images at a resolution of m were downsampled

by a factor of eight to arrive at images of size pixels and m

resolution Then for each potential mass object an ROI was extracted from

the corresponding downsampled image using the bounding box of the object

to dene the region A set of texture features based upon the GCMs yielded

a truepositive detection rate of at false positives per image and

detection accuracy at false positives per image with a dataset of

cases

Kobatake et al applied the iris lter for the detection of approximately

rounded convex regions and computed texture features based upon GCMs of

the iris lters output to isolate malignant tumors from normal tissue The

methods resulted in a detection rate of with false positives per

image with a dataset of CR images containing malignant tumors

Kegelmeyer developed a method to detect stellate lesions in mam

mograms and computed Laws texture features from a map of local edge

orientations A binary decision tree was used to classify the features Detec

tion results with ve test images yielded a sensitivity of with false

ndings per image Gupta and Undrill used Laws texture measures

and developed a texturebased segmentation approach to identify malignant

tumors

Researchers have also used fractal measures for the detection and char

acterization of mammographic lesions Priebe et al used texture and

fractal features for the detection of developing abnormalities Burdett et

al used fractal measures and nonlinear lters for the characterization of

lesion diusion Byng et al proposed measures of skewness of the image

brightness histogram and the fractal dimension of image texture which were

found to be strongly correlated with a radiologists subjective classications

of mammographic patterns

Gradientbased analysis of mass regions Several published algorithms

for breast mass segmentation are based upon the analysis of the gradient

orientation at edges to locate radiating lines or spicules

Karssemeijer and te Brake reported on the use of scale

space operators and linebased pixel orientation maps to detect spiculated

distortions The map of pixel orientations in the image was used to construct

two operators sensitive to starlike patterns of lines A total of images

nine spiculated masses cases of architectural distortion and normal

cases from the MIAS database were analyzed By combining the output

from both the operators in a classier an accuracy of was achieved in

detecting stellate carcinomas at one false positive per image The mass

cases studied belonged to the spiculated malignant category only and did not

include circumscribed malignant or circumscribed benign cases

In a related study te Brake and Karssemeijer used three pixel

based massdetection methods including the method they developed for de

tecting stellate distortions to examine if the detection of masses could be

achieved at a single scale Experiments with simulated masses indicated that

little could be gained by applying their methods at a number of scales Us

ing a dataset of images malignant and normals from the MIAS

database they reported that their gradient orientation method

performed better than two other methods based upon template

matching and Laplacian ltering Karssemeijer and te Brake mentioned

that their gradient orientation method requires the optimization of several

parameters in order to achieve a good detection accuracy

Kobatake and Yoshinaga developed skeleton analysis methods using

the iris lter to detect spicules of lesions They used a modied Hough trans

form to extract radiating lines from the center of a mass region to discrim

inate between starshaped malignant tumors and nonmalignant masses and

obtained an accuracy of in detecting malignant tumors using a dataset

of CR images including malignant nine benign and normal cases

Polakowski et al developed a modelbased vision algorithm using DoG

lters to detect masses and computed nine features based upon size circu

larity contrast and Laws texture features A multilayer perceptron neural

network was used for the classication of breast masses as benign or malig

nant With a dataset of malignant and benign cases they reported a

detection sensitivity of in identifying malignant masses with false

positives per image

Directional analysis of broglandular tissues Several researchers

have developed methods to analyze the ori

entation of broglandular tissues in order to detect mass regions and elim

inate false positives Zhang et al employed a Houghspectrumbased

technique to analyze the texture primitives of mammographic parenchymal

patterns to detect spiculated mass regions and regions possessing architec

tural distortion With a dataset of images obtained from patients the

methods yielded a sensitivity of with false positives per image Parr

et al studied Gabor enhancement of radiating spicules and concluded

that the linearity length and width parameters are signicantly dierent for

spicules in comparison to other linear structures in mammograms Later on

the group reported on the detection of linear structures in digital mammo

grams by employing a multiscale directional line operator Two

images indicating the probability of suspicion were derived by applying PCA

and factor analysis techniques to model the orientations present in central

core regions and the surrounding patterns of lesions respectively A sensi

tivity of at false positives per image was reported by combining

the evidence present in both the probability images in a knearestneighbor

approach using a dataset of mammograms containing spiculated lesions

and normal mammograms The basis for the computation of the

result mentioned above was not claried a value of false positives per

image with a dataset of mammograms leads to a total of less than one false

positive detection The sizes of the abnormalities in their studies ranged from

mm to mm with a mean of mm of the lesions tested were

smaller than mm Mudigonda et al applied techniques of texture

oweld analysis to analyze oriented patterns in mammograms by

computing the angle of anisotropy or dominant orientation and the strength

of owlike information or coherence at every point in the image the related

methods and results are described in the following subsections

Analysis of bilateral asymmetry Several studies

have investigated the bilateral asymmetry between the left and right

mammograms of an individual in order to localize breast lesions Such stud

ies have reported on the registration and unwarping transformations that are

required for the comparison of dierent images Apart from the anatomical

dierences that exist between the left and right breasts of an individual reg

istration methods will have to cope with the additional complexities arising

due to the variation in the amounts of compression applied in the process of

obtaining the images as well as the dierences in the angles of the projections

The transformations used in registration methods often face limitations in ad

equately addressing the abovementioned complexities of the mammographic

imaging process

Lau and Bischof applied a transformation based upon the outline of the

breast to identify asymmetry in breast architecture Using a Bspline model

of the breast outline to normalize images they compared features including

brightness roughness and directionality to dene asymmetry measures for

breast tumor detection Although they reported success in the detection of

tumors they warned that the method by itself is not reliable for clinical

application due mainly to the high rate of false negatives Giger et al

and Yin et al used similar warping procedures to align bilateral images

to perform subtraction Their method exploits the normal bilateral symmetry

of healthy parenchyma to label regions of signicant dierence as potential

masses Nishikawa et al analyzed asymmetric densities to detect masses

by performing nonlinear subtraction between the right and left breasts their

computeraided scheme for the classication of masses using articial neural

networks achieved a higher accuracy than radiologists

Miller and Astley proposed a method for the detection of asym

metry by comparing anatomically similar and homogeneous regions the pro

cedure included segmentation classication as fat or nonfat region texture

energy and shape features such as compactness circularity and eccentrici

ty Training and assessment were carried out on a leaveoneout basis with a

dataset of mammogram pairs achieving correct classication with a

linear discriminant classier Ferrari et al developed methods to achieve seg

mentation of the outline of the breast the pectoral muscle and the

broglandular disc in mammograms and analyzed the glandular tissue

patterns by applying directional ltering concepts to detect bilateral asym

metry Their methods described in Section avoid the application of

registration procedures that are associated with most of the abovementioned

methods to analyze asymmetry Instead they performed a global analysis of

the directional distributions of the glandular tissues using Gabor wavelets to

characterize the asymmetry in tissue ow patterns

Analysis of prior mammograms It is known that a small but sig

nicant number of cancer cases detected in screening programs have prompts

visible in earlier screening examinations These cases represent screening

errors or limitations and may occur due to the lack of an adequate under

standing of the perceptual features of early breast abnormalities as apparent

on mammograms The above observations have prompted many researchers

to analyze the previous or prior mam

mograms taken as part of routine screening and followup studies of patients

diagnosed with cancer in an eort to detect the disease in the prior mam

mograms Brzakovic et al developed methods for detecting changes in

mammograms of the veried positive group in the regions labeled as abnor

mal by a medical expert by comparing them with the features derived from

the corresponding regions of the previous screenings Their procedures in

cluded segmentation partitioning into statistically homogeneous regions and

regionbased statistical analysis in order to achieve registration and perform

comparison between mammograms acquired at dierent instances of screen

ing

Sameti et al studied the structural dierences between the re

gions that subsequently formed malignant masses on mammograms and other

normal areas in images taken in the last screening instance prior to the detec

tion of tumors Manually identied circular ROIs were transformed into their

optical density equivalents and further divided into three discrete regions

representing low medium and high optical density Based upon the discrete

regions a set of photometric and texture features was extracted They re

ported that in of the breast cancer cases studied it was possible to

realize the dierences between malignant mass regions and normal tissues in

previous screening images

Petrick et al studied the eectiveness of their massdetection

method in the detection of masses in prior mammograms The dataset used

included images malignant and benign from cases malignant

and benign Their detection methods achieved a by lm massdetection

sensitivity of with false positives per image They achieved a slightly

better accuracy of in detecting only malignant tumors Their detection

scheme attempts to segment salient densities by employing region growing

after enhancement of contrast in the image Such an intensitybased seg

mentation approach fails to detect the developing densities in the previous

screening images due to the inadequate contrast of mass regions before the

masses are actually formed

Several semiautomated segmentation schemes have

used manually segmented ROIs in order to search for masses in a specied

region of the breast Such methods nd limited practical utility in a screening

program

The female breast is a complex organ made up of brous glandular fatty

and lymphatic tissues The dierences in density information of the breast

tissues are captured in a mammogram in the form of intensity and textural

variations Mudigonda et al proposed an unsupervised segmentation

approach to localize suspicious mass regions in mammographic images The

approach aims to isolate the spatially interconnected structures in the image

to form regions concentrated around prominent intensities It would then be

possible to extract highlevel information characterizing the physical proper

ties of mass regions and to shortlist suspicious ROIs for further analysis A

block diagram of this approach is shown in Figure the various steps of

the detection algorithm are explained in detail in the following subsections

FIGURE

Block diagram of the massdetection algorithm Figure courtesy of NR

Mudigonda

Framework for pyramidal decomposition

Malignant tumors due to their invasive nature possess heterogeneous den

sity distributions and margins causing distortion in the orientation of the

surrounding tissues In order to detect such structures as single entities prior

smoothing of the image is required Mudigonda et al employed recursive

wavelet decomposition and Gaussian smoothing operations in a multiresolu

tion pyramidal architecture as preprocessing steps to achieve the required level

of smoothing of the image

A pyramidal representation of the given image was obtained by iterative

decimation operations on the fullresolution image thereby generating a hi

erarchy of subimages with progressively decreasing bandwidth and increasing

scale Wavelet decomposition divides the frequency spectrum of the

original image f into its lowpasssubbandequivalent image f

L

and highpass

equivalent detail image f

H

at dierent scales The lowpasssubband image

at each scale produced by decimating its preceding higherresolution image

present in the hierarchy by an octave level was further smoothed by a

Gaussian kernel and the resulting image was stretched to the range of

in pixel value The wavelet used was a symlet of eighth order Symlets are

compactly supported wavelets with the least asymmetry and the highest num

ber of vanishing moments for a given support width Figure shows

plots of the decomposition lowpass kernels used with symlets at two dierent

scales The wavelet decomposition was performed recursively to three octave

levels using the symlets mentioned above

The preprocessing steps of wavelet decomposition and Gaussian smooth

ing operations described above successively and cumulatively modulate the

intensity patterns of mass regions to form smooth hills with respect to their

surroundings in lowresolution images Figure a shows a

section of a mammogram containing two circumscribed benign masses Parts

b d of the gure show the corresponding lowresolution images after the

rst second and third levels of decomposition respectively The eects of the

preprocessing steps mentioned above may be observed in the lowresolution

images

The choice of the wavelet the width of the kernel used for lowpass ltering

and the degree or scale factor of decomposition can inuence the smoothed

results The preprocessing operations described above were employed to ar

rive at an estimate of the extent of isolated regions in a lowresolution image

and studies were not performed with dierent sets of choices However sat

isfactory smoothed results were obtained with the wavelet chosen due to its

symmetry A scale factor of three which causes the decomposition of the

original mpixel images to a resolution of mpixel was found to be

eective on most of the images tested decomposition to a higher scale resulted

in oversmoothing of images and merged multiple adjoining regions into single

large regions whereas a scale factor of two yielded insignicant regions due to

FIGURE

Plots of symlet decomposition lowpass lters at two scales Figure courtesy

of NR Mudigonda

a

b

c

d

FIGURE

a A section of a mammogram containing two circumscribed

benign masses Pixel size m Image width mm Lowresolution

images obtained by wavelet ltering b After the rst level of decomposition

pixels m per pixel c After two levels of decomposition

pixels m per pixel d After three levels of decomposition

pixels m per pixel The intensity of the ltered images has

been enhanced by four times for display purposes Figure courtesy of NR

Mudigonda

insucient smoothing However some researchers have performed mass

detection after reducing images to a resolution of mpixel

Segmentation based upon density slicing

The recursive smoothing and decimation operations described above result

in a gradual modulation of intensity information about the local intensity

maxima present in various isolated regions in the lowresolution image As a

result the intensity levels are expected to assume either unimodal or bimodal

histogram distributions The next step in the algorithm is to threshold the

image at varying levels of intensity to generate a map of isointensity con

tours The purpose of this step is to extract concentric groups of closed

contours to represent the isolated regions in the image

The densityslicing or intensityslicing technique slices the given image rep

resented as a D intensity function by using a plane that is placed parallel

to the coordinate plane of the image A level curve also known as

an isointensity curve is then formed by extracting the boundary of the area

of intersection of the plane and the intensity function Figure shows a

schematic illustration of the densityslicing operation Each level curve ob

tained using the procedure explained above is guaranteed to be continuous

and closed The number of levels of thresholding starting with the maxi

mum intensity in the image and the stepsize decrement for successive levels

were adaptively computed based upon the histogram distribution of the image

under consideration as explained below

Let f

max

represent the maximum intensity level in the lowresolution image

which was scaled to and let f

th

be the threshold representing the mass

tobackground separation which is to be derived from the histogram It

is assumed that the application of the preprocessing smoothing operations

results in exponentially decreasing intensity from the central core region of a

mass to its background represented as f

th

f

max

exp N where N is the

number of steps required for the exponentially decreasing intensity function to

attain the background level represented by f

th

N f

max

f

th

and is the

intended variation in step size between the successive levels of thresholding

The step size may be computed through a knowledge of the parameters f

th

and N The threshold f

th

was derived from the histogram and corresponds

to the intensity level representing the maximum number of occurrences when

the histogram assumes a unimodal distribution

It is essential to set bounds for f

th

so as not to miss the detection of masses

with lowdensity core regions while maintaining the computational time of

the algorithm at a reasonable level A large threshold value might miss the

grouping of lowintensity mass regions thereby aecting the sensitivity of the

detection procedure on the other hand a low value would result in a large map

of isointensity contours and increase the computational load on the algorithm

in further processing of the contours A smaller threshold could also result

in large numbers of false detections Initial estimates of f

th

derived from the

FIGURE

Schematic illustration of the densityslicing operation f

max

represents the

maximum intensity in the image and levels N represent a set

of N threshold values used for density slicing Figure courtesy of NR

Mudigonda

corresponding histograms of lowresolution images were observed to range

between and of f

max

and N was observed to range between and

The parameter f

th

was adaptively selected based upon the histogram as

explained below

If f

max

f

th

f

max

f

th

could be assumed to represent the mass

tobackground transition and the same threshold value is retained

If f

th

f

max

the mass regions that are to be detected in the image

are expected to be merged with the surrounding background and no

distinct central core regions would be present In such cases f

th

is

considered to be f

max

and N is set to the maximum number

of levels of thresholding considered to limit the stepsize increments of

the level function to a low value These steps facilitate close tracking of

diculttodetect masstobackground demarcation

If f

th

f

max

f

th

might not represent the true masstobackground

transition and hence is ignored An alternative search for f

th

is initi

ated so that the value obtained will lie in the upper half of the histogram

distribution

The steps described above realize a domain of isointensity contours in the

lowresolution image

Hierarchical grouping of isointensity contours

The next step in the algorithm is to perform grouping and elimination op

erations on the framework of closed contours generated in the lowresolution

image considering their parentchild nodal relations in a familytree architec

ture A schematic representation of such a hierarchical grouping procedure

is shown in Figure which depicts the segmentation of the lowresolution

image into three isolated regions based upon three concentric groups of con

tours

The strategy adopted was to shortlist at rst the possible central densecore

portions which are usually small in size but of higher density represented by

f

max

in each group of contours in Figure and to identify the immediate

lowdensity parent members encircling them The process was continued until

all the members in the available set of closed contours in the image were vis

ited Each of the closed contours was assigned to a specic group or family of

concentric contours based upon nodal relations thus leading to segmentation

of the image into isolated regions A concentric group of contours represents

the propagation of density information from the central core portion of an ob

ject in the image into the surrounding tissues In some images with dense and

fatty backgrounds the outermost contour members were observed to contain

multiple regions of dissimilar structures For this reason a specied number of

outer contours were discarded to separate the groups of contours representing

adjacent structures

The outermost contour in each family or group and the family count in

terms of the number of contours present could be useful in the analysis of

the regions segmented in order to reject false positives Masses irrespective

of their size were observed to result in a higher family count as compared

to elongated glandular tissues By setting a threshold on the family count

chosen to be ve dense glandular structures could be avoided from further

analysis A lower threshold value for the minimumallowable family count

was observed to aect the specicity in terms of the number of false positives

detected but not aect the sensitivity of the detection procedure Finally

the outermost contour from each of the shortlisted groups was upsampled to

the fullresolution image to form the corresponding segmented area

Results of segmentation of masses

The results of application of the algorithm to the image shown in Figure

a are presented in Figure a The contour map and the outer

most contours detected are shown superimposed on the lowresolution image

at a scale of three Figure shows the histogram of the correspond

ing lowresolution and smoothed image Figure b shows the contours

white upsampled from the lowresolution image in Figure a to the full

resolution image of Figure a the corresponding contours black that

were manually drawn by an expert radiologist are overlaid for comparison

FIGURE

Schematic representation of hierarchical grouping of contours G G and

G are groups of contours that represent isolated regions in the image Repro

duced with permission from NR Mudigonda RM Rangayyan and JEL

Desautels Detection of breast masses in mammograms by density slicing and

texture oweld analysis IEEE Transactions on Medical Imaging

c

IEEE

Figure shows a similar set of results of application of the massdetection

algorithm to a section of a mammogram containing a spiculated

malignant tumor

As can be seen from Figures and the upsampled contours in the

fullresolution image contain the most signicant portions of the correspond

ing masses and are in close agreement with the corresponding areas manually

delineated by the radiologist The results of application of the methods dis

cussed above to fullsize mammograms are presented in the subsections to fol

low along with methods based upon texture oweld principles for detailed

analysis of the various regions segmented in order to reject false positives

The massdetection algorithm was tested on segments of size up to

pixels of mammographic images benign and malignant from

the MIAS database with a spatial resolution of m m In

of the cases benign and malignant the segmented regions were

in agreement with the corresponding regions that were manually identied

by the radiologist In six images including ve images with circumscribed

benign masses and an image with a spiculated malignant tumor the regions

segmented by the algorithm were not in agreement with the corresponding

regions manually delineated by the radiologist The radiologist indicated that

he encountered diculty while tracing the boundaries of the masses in some

images from the MIAS database In the remaining four images where the

method failed all belonging to the spiculated benign category the mass por

tions are merged in fatty and glandular background In these images the

technique failed to generate a contour map with the specied minimum num

ber of concentric closed contours to be able to delineate the mass regions In

two of the images including a circumscribed benign mass and a spiculated

benign mass the masses are located close to the edges of the images and the

process of generation of concentric closed contours was impeded

Overall the massdetection algorithm performed well on images contain

ing malignant tumors and successfully segmented tumor areas that were in

agreement with the corresponding regions identied manually by the radi

ologist However the method encountered limited success in images with

benign masses In the detection scheme only the contour map generated in

the lowestresolution image is analyzed to segment the mass regions and the

information available in the intermediateresolution images along the hierar

chy is not considered Establishment of reliable intensity links through the

intermediateresolution images may result in improved detection results

Benignversusmalignant pattern classication was carried out using the

BMDP M stepwise discriminant analysis program with texture features

computed based upon averaged GCMs for the masses benign and

malignant that were successfully segmented by the massdetection procedure

See Sections and for details on the computation of texture features

using adaptive ribbons Four eective features including entropy second

moment second dierence moment and correlation were shortlisted The

a

Figure b

FIGURE

a Groups of isointensity contours and the outermost contour in each group in

the third lowresolution image of the mammogram section of Figure d

b The contours white of two masses indicated by arrows and two false

positives detected in the fullresolution image of Figure a with the

corresponding contours black of the masses drawn independently by a ra

diologist Reproduced with permission from NR Mudigonda RM Ran

gayyan and JEL Desautels Segmentation and classication of mammo

graphic masses Proceedings of SPIE Volume Medical Imaging

Image Processing pp

c

SPIE

FIGURE

Histogram of the lowresolution and smoothed image shown in Figure a

Reproduced with permission from NR Mudigonda RM Rangayyan and

JEL Desautels Segmentation and classication of mammographic masses

Proceedings of SPIE Volume Medical Imaging Image Processing

pp

c

SPIE

a

Figure b

c

Figure d

FIGURE

a A section of a mammogram containing a spiculated ma

lignant tumor Pixel size m Image width mm b Group of

isointensity contours and the outermost contour in the group in the third low

resolution image c Histogram of the lowresolution and smoothed image

shown d The contour white of the spiculated malignant tumor detected

in the fullresolution image superimposed with the corresponding contour

black drawn independently by a radiologist Reproduced with permission

from N R Mudigonda R M Rangayyan and J E L Desautels Segmenta

tion and classication of mammographic masses Proceedings of SPIE Volume

Medical Imaging Image Processing pp

c

SPIE

GCMbased texture features computed from the mass ribbons resulted in an

average classication eciency of

Detection of masses in full mammograms

Masses containing important signs of breast cancer may be dicult to detect

as they often occur in dense glandular tissue Successful identication of such

diculttodetect masses often results in a large number of false positives

Rejection of false positives forms an important part of algorithms for mass

detection

In the algorithm proposed by Mudigonda et al to detect masses the

pyramidal decomposition approach described in the preceding subsections was

extended for application to full mammograms Furthermore the orientation

information present in the margins of the regions detected was analyzed using

texture oweld principles to reject false positives The approach described

below is signicant in the following aspects

The methods constitute a comprehensive automated scheme for the de

tection of masses analysis of false positives and classication of mam

mographic masses as benign or malignant The detection methods pro

posed are not limitied to the analysis of any specic category of masses

instead they cover a wide spectrum of masses that include malignant tu

mors as well as benign masses of both the circumscribed and spiculated

categories

As described at the beginning of this section many of the recently

published research works have noted the

signicance of analyzing the oriented information in mammograms in

identifying regions that correspond to abnormal distortions in the im

ages Such methods require precise computation of orientation esti

mates Mudigonda et al introduced methods to analyze oriented

textural information in mammograms using oweld principles and

proposed features to dierentiate mass regions from other dense tissues

in the images in order to reduce the number of false positives The ow

eld methods employed use a strong analytical basis in order to provide

optimal orientation estimates

As discussed in Section the analysis of textural information in

ribbons of pixels across the boundaries of masses has been found to be

eective in the benignversusmalignant discrimination of masses

The studies cited above computed textural features using

ribbons of pixels of a xed width of mm In the work of Mudigonda et

al a method is proposed to estimate the widths of the ribbons to

be able to adapt to variations in the size and shape of the detected mass

regions The adaptive ribbons of pixels extracted are used to classify the

regions detected as masses or false positives at rst and subsequently

to discriminate between benign masses and malignant tumors

The features used for the classication of masses and false positives are

based upon specic radiographic characteristics of masses as described

by an expert radiologist

The block diagram shown in Figure lists the various steps of the de

tection algorithm which are explained in detail in the following paragraphs

Detection of the breast boundary In order to limit further processing

to only the breast region an approximate outline of the breast was detected

initially by employing the following steps The image was smoothed with a

separable Gaussian kernel of width pixels pixel width m and

quantized to gray levels The method proposed by Schunck and Rao

and Schunck was used to generate Gaussian kernels Figure shows

a plot of the Gaussian kernel used

A map of isointensity contours was generated by thresholding the image

using a threshold close to zero From the map of isointensity contours a

set of closed contours was identied by employing the chain code The

contour containing the largest area was then considered to be the outline of

the breast Figure illustrates a mammogram of size pixels

with a spiculated malignant tumor The outline of the breast detected in the

mammogram of Figure is shown in Figure The method successfully

detected the outlines of all of the images tested The method worked

successfully with images lacking skinair boundaries all around as well

Detection of salient densities Gaussian pyramidal decomposition was

employed to achieve the required smoothing instead of wavelet decomposition

that was used in the case of detection of masses in sectional images of mam

mograms as discussed in Section Gaussian decomposition when tested

with some of the previously used sectional mammographic images provided

comparable detection results in terms of the boundaries of the regions seg

mented Gaussian smoothing using separable kernels has the advantage of

ease of implementation and also provides computational advantages partic

ularly with fullsize mammograms

FIGURE

Block diagram of the algorithm for the detection of masses in full mammo

grams Reproduced with permission from N R Mudigonda R M Rangayyan

and J E L Desautels Detection of breast masses in mammograms by den

sity slicing and texture oweld analysis IEEE Transactions on Medical

Imaging

c

IEEE

FIGURE

Plot of a Gaussian kernel with the support width of pixels The width at

halfmaximum height is ve pixels Figure courtesy of N R Mudigonda

The original b images with a spatial resolution of m were subsam

pled to a resolution of m after performing smoothing with a separable

Gaussian kernel of width ve pixels The width of the Gaussian kernel at

halfmaximum height is about m and hence is not expected to cause

excessive smoothing of mass regions because mass features in mammograms

typically span a few millimeters

The preprocessing steps described above are essential in order to capture

the complete extent of mass features as single large regions so as to facilitate

adequate inputs for further discrimination between masses and false positives

Masses were assumed to be hyperdense or at least of the same density with

respect to their background The preprocessing steps described above may

not be eective with masses that do not satisfy this assumption

Multilevel thresholding In the procedure of Mudigonda et al

the lowresolution image is initially reduced to gray levels in intensity

and thresholded at N levels starting from the maximum intensity level

f

max

with a stepsize decrement of f

max

The purpose of this

step is to extract concentric groups of closed contours to represent the isolated

regions in the image as explained earlier The abovementioned parameters

were chosen based upon the observation of the histograms of several low

resolution images The histogram of the lowresolution image obtained by

preprocessing the mammogram in Figure is shown in Figure As

indicated in the gure the intensity level at which the masses and other dense

FIGURE

A mammogram size pixels m per pixel with a spic

ulated malignant tumor radius cm Case mdb from the MIAS

database Reproduced with permission from N R Mudigonda R M

Rangayyan and J E L Desautels Detection of breast masses in mammo

grams by density slicing and texture oweld analysis IEEE Transactions

on Medical Imaging

c

IEEE

FIGURE

The map of isointensity contours extracted in the smoothed and subsampled

version size pixels m per pixel of the mammogram shown in

Figure The breast outline detected is superimposed In some cases sev

eral contours overlap to produce thick contours in the printed version of the

image Reproduced with permission from N R Mudigonda R M Rangayyan

and J E L Desautels Detection of breast masses in mammograms by den

sity slicing and texture oweld analysis IEEE Transactions on Medical

Imaging

c

IEEE

tissues appear to merge with the surrounding breast parenchyma is around

the minimum threshold level of

FIGURE

Histogram of the lowresolution image corresponding to the mammogram in

Figure Reproduced with permission from N R Mudigonda R M Ran

gayyan and J E L Desautels Detection of breast masses in mammograms by

density slicing and texture oweld analysis IEEE Transactions on Medical

Imaging

c

IEEE

Figure shows the map of isointensity contours obtained by density

slicing or multilevel thresholding of the mammogram shown in Figure

Observe that multiple concentric contours may fuse into thick contours in

the printed version of the image

Grouping of isointensity contours The scheme represented in Fig

ure was adopted to perform a twostep grouping and merging operation

on the individual contours possessing a minimum circumference of mm ve

pixels at m to arrive at groups of concentric isointensity contours

Initially the contour members with intensity values ranging from f

max

to

f

max

with f

max

were grouped to form a set of regions corresponding

to high intensities in the image and then the remaining contour members

were grouped into a separate set The undesired merging of adjoining regions

was controlled by monitoring the running family count of each group for any

abrupt uctuations in terms of its family count The information from both

the sets of groups of contours was combined by establishing correspondences

among the outermost members of the various groups present in each set to

arrive at the nal set of segmented regions in the lowresolution image The

largest contour in each group thus nalized with a minimum family count

of two members was upsampled into the fullresolution image to form the

corresponding segmented area The segmented regions identied in the full

resolution image were analyzed using the features described above to identify

truepositive and falsepositive regions

Analysis of mammograms using texture oweld

In a mammogram of a normal breast the broglandular tissues present ori

ented and owlike or anisotropic textural information Mudigonda et al

proposed features to discriminate between masses and the strongly oriented

broglandular tissues based upon the analysis of oriented texture in mammo

grams The method proposed by Rao and Schunck briey described

in the following paragraphs was used to characterize owlike information in

the form of intrinsic orientation angle and coherence images The intrinsic

angle image reveals the direction of anisotropy or ow orientation of the un

derlying texture at every point in the image Coherence is a measure of the

degree or strength of anisotropy in the direction of ow

Rao and Schunck and Rao made a qualitative comparison of the

performance of their method with the method proposed by Kass and Witkin

and reported that their method achieved superior results in character

izing oweld information However it appears that their implementation

of Kass and Witkins scheme diered from the method that was originally

proposed by Kass and Witkin Regardless the method of Rao and Schunck

has a strong analytical basis with fewer assumptions

The methodology to derive the intrinsic images begins with the computation

of the gradient information at every point in the image by preprocessing the

image with a gradientofGaussian lter of a specied width The impulse

response of a D Gaussian smoothing lter gx y of width is given by

gx y exp

x

y

where the scale factor has been ignored The impulse response of the gradient

ofGaussian lter hx y tuned to a specied orientation can be obtained

using gx y as

hx y

g

x

g

y

cos sin

where represents the dot product At each point in the given image the

lter hx y upon convolution with the image yields the maximal response in

the orientation that is perpendicular to the orientation of the underlying

texture that is the angle of anisotropy Based upon the above and with

the assumption that there exists a dominant orientation at every point in the

given image Rao and Schunck derived the optimal solution to compute

the angle of anisotropy

pq

at a point p q in the image as described below

LetG

mn

and

mn

represent the gradient magnitude and gradient orientation

at the point mn in an image respectively and P P be the size of the

neighborhood around p q used for computing

pq

The gradient magnitude

is computed as

G

mn

q

G

x

mn G

y

mn

where G

x

mn and G

y

mn represent the outputs of the gradientofGauss

ian lter at mn in the x and y directions respectively The gradient ori

entation is computed as

mn

arctan

G

y

mn

G

x

mn

The projection of G

mn

on to the gradient orientation vector at p q at angle

pq

is G

mn

cos

mn

pq

as illustrated schematically in Figure

Based upon the discussion above the sumofsquares S of the projections of

the gradient magnitudes computed at the various points of the neighborhood

in a reference orientation specied by is given by

S

P

X

m

P

X

n

G

mn

cos

mn

The sum S varies as the orientation is varied and attains its maximal

value when is perpendicular to the dominant orientation that represents the

underlying texture in the given set of points Dierentiating S with respect

to yields

dS

d

P

X

m

P

X

n

G

mn

cos

mn

sin

mn

By setting

dS

d

and further simplifying the result we obtain the solution

for

pq

that maximizes S at the point p q in the image as

pq

arctan

P

P

m

P

P

n

G

mn

sin

mn

P

P

m

P

P

n

G

mn

cos

mn

The second derivative

d

S

d

is given by

d

S

d

P

X

m

P

X

n

G

mn

cos

mn

FIGURE

Schematic illustration of the projection of the gradient magnitude for com

puting the dominant orientation angle and coherence the scheme of Rao and

Schunck G

pq

and

pq

indicate the gradient magnitude and orientation

at p q respectively The corresponding parameters at mn are G

mn

and

mn

The size of the neighborhood shown is P P pixels Figure

courtesy of N R Mudigonda

The value of

pq

that is obtained using Equation represents the direc

tion of the maximal gradient output because the second derivative shown in

Equation is negative at

pq

when the texture has only one dominant

orientation The estimated orientation angle of ow

pq

at p q in the image

is then

pq

pq

because the gradient vector is perpendicular to the direction of ow The

angles computed as above range between and radians

In order to analyze mammograms with the procedure described above the

original image was initially smoothed using a separable Gaussian kernel

of a specied width and the gradients in the x and y directions were computed

from the smoothed image using nite dierences in the respective directions

The choice of the width of the Gaussian aects the gradient computation a

width of mm pixels was used by Mudigonda et al in relation to

the range of the size of features related to breast masses The lter has a width

of about mm at its halfmaximum height This lter size is appropriate

given that mammograms may demonstrate lumps that are as small as mm

in diameter

The gradient estimates computed as above were smoothed using a neigh

borhood of size pixels mm the width of which was chosen to

be larger than the Gaussian that was initially used to compute the gradient

estimates Figure shows the intrinsic angle image of the mammogram

shown in Figure The bright needles overlaid on the image indicate the

underlying dominant orientation at points spaced every fth row and fth

column computed using Equation Needles have been plotted only for

those pixels where the coherence computed as follows is greater than zero

The coherence

pq

at a point p q in the given image was computed as the

cumulative sum of the projections of the gradient magnitudes of the pixels in

a window of size P P in the direction of the dominant orientation at the

point p q under consideration as

pq

G

pq

P

P

m

P

P

n

G

mn

cos

mn

pq

P

P

m

P

P

n

G

mn

The result was normalized with the cumulative sum of the gradient magni

tudes in the window and multiplied with the gradient magnitude at the point

under consideration in order to obtain high coherence values at the points in

the image having high visual contrast The coherence image computed for

the mammogram in Figure is shown in Figure It can be observed

that glandular tissues ligaments ducts and spicules corresponding to archi

tectural distortion possess high coherence values

FIGURE

Intrinsic angle information white lines for the mammogram shown in Fig

ure The boundaries black represent the mass and falsepositive regions

segmented at the initial stage of the massdetection algorithm The breast

outline detected is superimposed Reproduced with permission from N R

Mudigonda R M Rangayyan and J E L Desautels Detection of breast

masses in mammograms by density slicing and texture oweld analysis

IEEE Transactions on Medical Imaging

c

IEEE

FIGURE

Intrinsic coherence image of the mammogram shown in Figure Repro

duced with permission from N R Mudigonda R M Rangayyan and J E L

Desautels Detection of breast masses in mammograms by density slicing and

texture oweld analysis IEEE Transactions on Medical Imaging

c

IEEE

Adaptive computation of features in ribbons

The regions detected by the method described above vary greatly in size and

shape For this reason a method was devised to compute adaptively the width

of the ribbon for the derivation of features see Section or equivalently

the diameter of the circular morphological operator for a particular region

based upon the regions size and shape

Figure shows a schematic representation of the method used to com

pute adaptively the size of the ribbon Initially the diameter of the bounding

circle enclosing a given candidate region was found by computing the maxi

mal distance between any two points on its boundary Then the areas of the

region A

r

and the bounding circle A

c

enclosing the region were computed

The width of the ribbon was computed as

R

w

R

c

A

r

A

c

where R

c

is the radius of the bounding circle The ratio

A

r

A

c

is a simple measure

of the narrowness and shape complexity of the region The size of the ribbon

computed above was limited to a maximum of mm or pixels The regions

for which the sizes of ribbons computed was less than mm or four pixels

were rejected and not processed further in the falsepositive analysis stage

The ribbons of pixels white extracted across the boundaries black of the

various regions detected in the image shown in Figure are illustrated in

Figure

FIGURE

Schematic representation of the adaptive computation of the width of the

ribbon A

r

area of the candidate region A

c

area of the bounding circle and

R

c

radius of the bounding circle Figure courtesy of N R Mudigonda

FIGURE

Ribbons of pixels white extracted adaptively across the boundaries black

of the regions detected in the mammogram shown in Figure Repro

duced with permission from N R Mudigonda R M Rangayyan and J E L

Desautels Detection of breast masses in mammograms by density slicing and

texture oweld analysis IEEE Transactions on Medical Imaging

c

IEEE

Features for massversusfalsepositive classi cation In order to

classify the regions detected as true masses or false positives the following

features were proposed by Mudigonda et al based upon certain well

established radiological notions about breast masses as apparent on mammo

grams

Contrast C

fg

Masses in mammograms may be presumed to be hy

perdense or at least isodense with respect to their surroundings For

this reason the contrast C

fg

of a region was computed as the dier

ence between the mean intensities of the foreground region or ROI and

a background region dened as the region enclosed by the extracted rib

bon of pixels excluding the ROI Regions possessing negative contrast

values were rejected from further analysis

Coherence ratio

r

The interior regions of masses are expected to be

less coherent than their edges The ratio

r

of the mean coherence of

the ROI excluding the ribbon of pixels to the mean coherence in the

ribbon of pixels was used as a feature in pattern classication

Entropy of orientation estimates H

o

The orientation of spicules in

the margins of spiculated masses is usually random Furthermore the

orientation estimates computed in the margins of circumscribed masses

could cover a wide range of angles between zero and radians and may

not possess any dominant orientation On the contrary broglandular

tissues are highly directional For these reasons the entropy H

o

of

the orientation estimates was computed in the ribbon of pixels of each

region detected for use as a feature in pattern classication

Variance of coherenceweighted angle estimates

h

The fourth fea

ture was based upon the coherenceweighted angular histogram which

was computed for a particular region by incrementing the numbers of

occurrence of angles with the magnitudes of coherence values computed

at the respective points after resampling the angle values in the rib

bon regions to Q equally spaced levels between zero and This

is equivalent to obtaining a cumulative sum of the coherence estimates

of the points belonging to each bin as the height of the corresponding

bin in the histogram The histogram distributions obtained as above

were normalized with the cumulative sum of the coherence values com

puted in the ribbons of the respective regions and the variance

h

was

computed as

h

Q

Q

X

i

i

h

where i i Q are the normalized values of the heights of

the Q bins of the histogram formed by the coherenceweighted angle

estimates and

h

is the average height of the bins of the histogram

h

Q

Q

X

i

i

Features for benignversusmalignant classi cation The ecacy in

benignversusmalignant classication of the truepositive mass regions suc

cessfully segmented by the massdetection algorithm was evaluated by using

a set of ve GCMbased texture features entropy second moment dier

ence moment inverse dierence moment and correlation see Section

for details The features were computed in the ribbon of pixels extracted

adaptively from each segmented mass margin as described above The GCMs

constructed by scanning each mass ribbon in the

o

o

o

and

o

di

rections were averaged to obtain a single GCM and the ve texture features

were computed for the averaged GCM for each ribbon

Results of mass detection in full mammograms

Mudigonda et al tested their methods with a total of images each of

size pixels at a resolution of m including benign masses

malignant tumors and normal cases selected from the MiniMIAS

database The dataset included circumscribed and spiculated cases in both

of the benign and malignant categories The mean values of the sizes of the

masses were cm and cm for the benign and malignant

categories respectively The radius of the smallest mass malignant was

cm and that of the largest mass benign was cm The center of

abnormality and an approximate radius of each mass are indicated in the

database The circular demarcation of masses as done in the database is not

useful for conrming the results of mass detection because such a demarcation

may also include normal broglandular tissues in the ROIs particularly in

spiculated cases Hence only the center of the abnormality as indicated for

each mass in the database was used to conrm the result of mass detection

The massdetection algorithm successfully detected all of the malignant

tumors in the database used However the algorithm met with limited success

in detecting benign masses In ve circumscribed and six spiculated of

the benign cases tested the algorithm failed to detect the masses The

overall detection accuracy was with a total of cases

A close inspection of some of the benign cases in which the algorithm failed

to detect the mass revealed the following details In three of the four dense

glandular masses cases labeled as mdb mdb and mdb in the Mini

MIAS database two fattyglandular masses mdb and mdb

and a fatty mass mdb the masses do not have prominent central core

regions and possess poor contrast with respect to their background Con

trast enhancement of such mammograms prior to the detection step

could improve the performance of the detection algorithm In a mammo

gram containing a fatty mass mdb the high intensity in the pectoral

muscle region aected the multilevel thresholding process of the detection al

gorithm This calls for methods to detect precisely the pectoral muscle

see Section and a twostep detection procedure initially the mammo

graphically apparent regions corresponding to suspicious lymph nodes could

be searched for inside the pectoral muscle region and in the second stage the

region of the breast excluding the pectoral muscle area could be searched for

the possible presence of masses

In two other cases mdb and mdb the masses did not satisfy the

hyperdense or isodense assumption Successful detection of masses in such

cases may require additional methods based upon the asymmetry between

the right and the left mammograms in order to detect regions possessing

architectural distortion

In the method proposed by Mudigonda et al no region was rejected

based upon its size during the initial stage because masses can be present

with any size Instead the emphasis was on reducing false positives by using

texture oweld features As a result a large number of falsepositive regions

at approximately per image were detected by the algorithm along with the

true mass regions during the initial stage of detection

Massversusfalsepositive classi cation The four features C

fg

r

H

o

and

h

described in Section were computed in the ribbons of the

candidate regions that were detected in all of the cases tested and used in

a linear discriminant classier to identify the true mass regions and false pos

itives The MIAS database contains an unusual number of spiculated benign

cases In order to study the eect of such an atypical distribution of cases

on the accuracy of the detection method pattern classication experiments

were carried out in two stages at rst a massversusnormaltissue classica

tion was conducted with the regions detected in the cases tested Next

malignanttumorversusnormaltissue classication was performed using the

features computed from the regions detected in the malignant and the

normal cases tested

Pattern classication was carried out using the BMDP M stepwise dis

criminant analysis program with the leaveoneout scheme In datasets

with limited numbers of cases it is dicult to form separate sets of data

for training and testing purposes The leaveoneout crossvalidation scheme

helps to obtain the leastbiased estimates of classication accuracy in such

situations The overall classication eciency in the classication of malig

nant tumors versus normal tissue was and that for discriminating between

masses both benign and malignant and normal tissue was

The linear discriminant function obtained at equal prior probability values

for the group with tumors and the group with normal cases was encoded

to arrive at a massversusfalsepositive decision for each segmented region

Figure indicates the nal set of regions that were retained for the image in

Figure after the detection and falsepositive analysis stages The region

detected inside the pectoral muscle area as shown in Figure was suspected

to be a lymph node aected by the invasive carcinoma Such areas may be

prompted to a radiologist for studying possible nodal involvement particularly

in cases with no localizing signs of the disease

Figure shows a mammogram with a spiculated malignant tumor in

dicated by an arrow radius cm that is smaller and less obvious than

the tumor shown in Figure Figures and show the results of de

tection for the image shown in Figure before and after the falsepositive

analysis stage respectively The tumor has been successfully detected along

with one false positive

The massversusnormaltissue classication experiment involving the

mass regions benign and malignant that the algorithm successfully

detected and false positives from a total of images including normal

cases resulted in an overall classication eciency of with a sensitivity

of at false positives per image A total of six masses four benign and

two malignant were misclassied as normal tissue However if the fact that

the algorithm missed benign masses during the initial stage of detection

itself is taken into consideration the true detection sensitivity of the algorithm

with the database of benign and malignant masses reduces to

In the case of malignanttumorversusnormaltissue classication a high

overall classication eciency of was achieved the dataset included

malignant tumors and false positives from a total of images including

normal cases A sensitivity of was obtained at false positives per

image Although all of the tumors were successfully detected in the initial

stage two of the malignant tumors that were detected were misclassied later

as normal tissue yielding a small proportion of false negatives

In a related work te Brake and Karssemeijer compared three mass

detection schemes to detect malignant tumors of small size radius smaller

than mm medium size radius between mm and mm and large size

radius greater than mm The method based upon gradientorientation

maps was reported to have achieved the best results in the detection

of masses in the small and medium categories The study of te Brake and

Karssemeijer focused only on malignant tumors and did not include benign

masses

Petrick et al reported similar trends in detection results using a

dataset of mammograms containing benign and malignant cases

A sensitivity of was reported at false positives per image Most

of the other previous studies in the related eld reported on

the detection of only malignant tumors and did not specically consider the

detection of benign masses

Benignversusmalignant classi cation The eectiveness of the seg

mentation results in benignversusmalignant pattern classication was veri

ed using the ve GCMbased texture features computed based upon aver

aged GCMs as explained above for the cases benign and malignant

FIGURE

Adaptive ribbons of pixels white and boundaries black of the regions re

tained in the mammogram shown in Figure after the falsepositive anal

ysis stage The larger region corresponds to the malignant tumor the other

region is a false positive See also Figure Reproduced with permission

from N R Mudigonda R M Rangayyan and J E L Desautels Detection

of breast masses in mammograms by density slicing and texture oweld

analysis IEEE Transactions on Medical Imaging

c

IEEE

FIGURE

A mammogram size pixels m per pixel with a spiculated

malignant tumor pointed by the arrow radius cm Case mdb from

the MIAS database Reproduced with permission from N R Mudigonda

R M Rangayyan and J E L Desautels Detection of breast masses in mam

mograms by density slicing and texture oweld analysis IEEE Transac

tions on Medical Imaging

c

IEEE

FIGURE

Ribbons of pixels white extracted adaptively across the boundaries black

of the regions detected in the mammogram shown in Figure Repro

duced with permission from N R Mudigonda R M Rangayyan and J E L

Desautels Detection of breast masses in mammograms by density slicing and

texture oweld analysis IEEE Transactions on Medical Imaging

c

IEEE

FIGURE

Adaptive ribbons of pixels white and boundaries black of the regions re

tained in the mammogram shown in Figure after the falsepositive anal

ysis stage The larger region corresponds to the malignant tumor the other

region is a false positive See also Figure Reproduced with permission

from N R Mudigonda R M Rangayyan and J E L Desautels Detection

of breast masses in mammograms by density slicing and texture oweld

analysis IEEE Transactions on Medical Imaging

c

IEEE

that were successfully segmented by the massdetection procedure Pattern

classication was carried out using the BMDP stepwise logistic regression

program The ve GCMbased texture features resulted in an overall

classication eciency of The results obtained conrm that the mass

regions segmented in images of resolution m possess adequate discrim

inant information to permit their classication as benign or malignant with

texture features Similar benignversusmalignant classication results were

obtained using partial images of the same cases but with m resolution

It appears that the detection and classication of masses may be successfully

performed using images of resolution m with the techniques described

above

Application Bilateral Asymmetry in Mammograms

Asymmetry between the left and right mammograms of a given subject is an

important sign used by radiologists to diagnose breast cancer Analysis

of asymmetry can provide clues about the presence of early signs of tumors

parenchymal distortion small asymmetric bright spots and contrast etc

that are not evaluated by other methods Several works have been

presented in the literature addressing this problem

with most of them applying some type of alignment of the breast images before

performing asymmetry analysis However alignment procedures applied to

mammograms have to confront many dicult problems such as the natural

asymmetry of the breasts of a given subject the absence of good corresponding

points between the left and right breast images to perform matching and the

distortions inherent to breast imaging

Procedures for systematic analysis were proposed by Lau and Bischof

and Miller and Astley to perform comparison of the corresponding

anatomical regions between the left and right breast images of an individual

in terms of shape texture and density Lau and Bischof also proposed

a directional feature to quantify oriented patterns

Ferrari et al proposed a procedure based upon directional analysis

using Gabor wavelets in order to analyze the possible presence of global dis

turbance between the left and right mammograms of an individual in the

normally symmetrical ow of mammary structures The analysis was focused

on the broglandular disc of the mammograms segmented in a preprocessing

step The methods and results of Ferrari et al are presented in the

following paragraphs

The broglandular disc

As indicated by the proceedings of the recent International Workshops on

Digital Mammography several researchers are developing

image processing methods to detect early breast cancer Most of the tech

niques proposed perform analysis of the whole mammogram without taking

into account the fact that mammograms have dierent density patterns and

anatomical regions that are used by radiologists in diagnostic interpretation

In fact mammographic images are complex and dicult to analyze due to

the wide variation in the density and the variable proportion of fatty and

broglandular tissues in the breast Based upon these observations a

few researchers have proposed methods to segment and also to model mam

mograms in terms of anatomical regions

Miller and Astley investigated the visual cues utilized by radiologists

and the importance of a comparison of the corresponding anatomical struc

tures in order to detect asymmetry between the left and right mammograms

of an individual However in their work the anatomical segmentation ap

proach and the possible methodologies to segment the mammograms were

not the main issue Aylward et al devised a modeling system to seg

ment a given mammographic image into ve major components background

uncompressed fat periphery of the breast close to the skinair boundary

fat dense tissue and muscle the system combined geometric and statistical

techniques A few other segmentation techniques for mammograms have been

presented in the literature however the focus has not been on anatomical seg

mentation but on specic problems such as density correction of peripheral

breast tissue localization of the nipple on mammograms

and quantication of breast density and its association with the risk of

breast cancer

The broglandular disc is an anatomical region of the breast characterized

by dense tissues ligaments and milk ducts Normally it has the shape of a

disc or a cone and goes through the interior of the breast from the region

near the chest wall to the nipple Segmentation of the broglandular

disc could form an important stage in techniques for the detection of breast

cancer that use asymmetry between the left and right mammograms of the

same subject or for monitoring breast density changes in screening programs

According to Caulkin et al it has been noticed clinically that breast

cancer occurs most frequently in the upper and outer quadrant of the breast

and that the majority of cancers are associated with glandular rather than

fatty tissues A common procedure used by radiologists in screening programs

for the detection of breast cancer is comparison between the left and right

broglandular discs of the mammograms of the same subject

Several works reported in the literature have been directed to address the

problem of automatic quantication of breast density and its association with

the risk of breast cancer Most of such

works propose an index or a set of values for the quantication of breast tissue

density However only a few works have attempted to address the problem of

detection and segmentation of the broglandular disc

for subsequent analysis

Ferrari et al proposed a method to segment the broglandular disc

in mammograms In this method prior to the detection of the broglan

dular disc the breast boundary and the pectoral muscle are detected using

the methods described in Sections and The broglandular disc is

detected by dening a breast density model The parameters of the model

are estimated by using the EM algorithm and the minimumdescription

length MDL principle Then a reference value computed by using

information from the pectoral muscle region is used along with the breast

density model in order to identify the broglandular disc The details of the

methods are described in the following subsections

Gaussian mixture model of breast density

The breast density model used by Ferrari et al is based upon a Gaussian

mixture model see also Section estimated by using the gray

level intensity distribution that represents categories or classes with dierent

density values in mammograms Except for the rst category the categories

are related to types of tissue that may be present in the breast

Dierent from the model proposed by Aylward et al which xes

at ve the number of tissue classes the model used by Ferrari et al

was formulated with the hypothesis that the number of tissue classes in the

eective region of the breast after extracting the pectoral muscle may vary

from two to four among the following possibilities

Uncompressed fatty tissues represented by fatty tissues localized in

the periphery of the breast

Fatty tissues composed by fatty tissues that are localized next to

the uncompressed fatty tissues and surround the denser areas of the

broglandular disc

Nonuniform density tissues including the density region that sur

rounds the highdensity portions of the broglandular disc extending

close to the chest wall

Highdensity tissues represented by the highdensity portions of the

broglandular disc

The hypothesis described above is based upon the fact that breast tissues

may naturally vary from one person to another or even for the same person

during her life time due to aging or hormonereplacement therapy

A highdensity and highfat breast for example will likely present only two

categories It was assumed that the data the graylevel values in a segmented

mammogram are generated by a Gaussian mixture model with a onetoone

correspondence between the mixture model components and the observed

tissue classes see also Section Thus the marginal probability of having

a gray level x is the sum of the probability over all mixture components and

it is represented by a linear superposition of multiple weighted Gaussians as

pxj

K

X

i

W

i

pxj

i

where the x values represent the graylevel values in the image W

i

are

the normalized mixing parameters

P

K

i

W

i

with W

i

pxj

i

is the Gaussian PDF parameterized by

i

i

i

that is the

mean value

i

and the standard deviation

i

of the i

th

Gaussian kernel the

vector represents the collection of the parameters of the mixture model

W

W

W

K

K

and K is the number of Gaussian kernels

that is tissue categories The Gaussian kernel is represented as

pxj

i

p

i

exp

x

i

i

In the case of using features other than the graylevel values of the image

such as texture features a multivariate Gaussian must be used instead of a

univariate Gaussian In this case the mean value and the standard deviation

of the graylevel values are replaced by the mean vector and the covariance

matrix of the feature vectors respectively In the model as above the Bayesian

assumption is made that the PDF associated with a pixel in the image is

independent of that of the other pixels given a class of tissue and furthermore

independent of its position in the image The estimation of the parameters is

performed by using the EM algorithm which is an iterative procedure that

maximizes the loglikelihood of the parameters of the model for a dataset

representing a PDF In the EM algorithm the estimation of the model

parameters is performed in two consecutive steps the Estep and the Mstep

In the Estep the current set of parameters is used to compute the model

The model is then assumed to be correct and the most likely distribution of

the data with respect to the model is found In the Mstep the parameters

of the model are reevaluated with respect to the new data distribution by

maximizing the loglikelihood given as

logLj log

N

Y

i

px

i

j

where N is the number of pixels in the eective region of the breast which is

the region demarcated by the breast boundary without the pectoral muscle

and represents the data sample The procedure is iterated until the values

of logLj between two consecutive estimation steps increase by less than

or the number of iterations reaches a specied limit cycles

Initialization of the model parameters The parameters of the model

were rst initialized by setting the center and weight of each Gaussian as

i

and W

i

K where i K is the index of the Gaussian

kernel and is a random value within the range dened by the minimum

and maximum graylevel values present in the eective area of the breast

The variance

i

of each Gaussian was initialized to the nearest distance to

the other Gaussian kernels If the variance

i

became less than unity during

the maximization step the Mstep it was reinitialized with a large random

value This procedure was intended to avoid shrinkage of the variance to a

small value The EM estimation procedure was initialized and repeated three

times in order to minimize the chance of convergence to a local minimum

Model selection Besides the initialization of the parameters another

diculty with the mixture model lies in choosing the number of components

that best suits the number of clusters or density categories present in the

image Several methods have been proposed in the literature to estimate the

number of components in a mixture of Gaussians among which Akaikes in

formation criterion Bayesian information criterion and MDL are the most

commonly used methods These estimators are motivated from dif

ferent theories that translate in practice to dierent penalty factors in the

formulation used to select the best model The MDL criterion is based upon an

informationtheoretic view of induction as data compression It is equivalent

to the Bayesian information criterion which gives a Bayesian interpretation

Akaikes information criterion is derived from a dierent theoretical perspec

tive it is an optimal selection rule in terms of prediction error that is the

criterion identies a nitedimensional model that while approximating the

data provided has good prediction properties However Akaikes information

criterion tends to yield models that are large and complex In general the

MDL criterion outperforms Akaikes and Bayesian criteria

As discussed above the number of density categories in the breast density

model can vary from two to four Because no reliable prior information is avail

able about the number of tissue categories present in a given mammogram

the MDL principle was used to select the number K of the Gaussian

kernels of the model The MDL principle deals with a tradeo between the

maximumlikelihood or minimumerror criterion for tting the model to a

dataset and the complexity of the model being designed Thus if the

models designed by using two dierent values of K t the data equally well

the simpler model is used The value of K was chosen so as to maximize the

quantity

logLj

NK

logK

where NK Kd is the number of free parameters in the mixture

model with K Gaussian kernels The value of K ranges from two to four and

d represents the dimension of the feature space

Delimitation of the broglandular disc

After computing the parameters of the Gaussian mixture model the maximum

likelihood method was applied to the original image to produce a Klevel

image that encoded at each pixel a cluster membership with the highest

likelihood among the K estimated Gaussian kernels Figure shows the

eective region of the image mdb used for the model estimation process

and the frequency distribution plots of the resulting Gaussian mixture com

ponents

According to Karssemeijer the density of the pectoral muscle is an

important item of information that can be used as a reference in the interpreta

tion of densities in the breast tissue area because regions of similar brightness

or density will most likely correspond to broglandular tissue Based upon

this observation and the breast density model described above a postprocess

ing stage was developed in order to determine the cluster region in the Klevel

image if one existed that agreed with the broglandular disc In this stage

the Klevel cluster was classied as the broglandular region if

K

P

P

where

P

and

P

are respectively the mean and standard deviation of the

graylevel values of the pectoral muscle region and

K

is the mean gray level

of the cluster K computed from the eective region of the given image The

threshold value

P

P

used to determine the broglandular disc was ar

rived at based upon experiments because a direct comparison between the

densities of the pectoral muscle and the broglandular tissue in terms of

physical parameters would be dicult due to several inuential factors of the

imageacquisition protocol

Figures e and f illustrate respectively the Klevel image K

by MDL resulting from the mixture model designed and the broglandular

disc identied for the mammogram in Figure a The results for another

mammogram are provided in Figure with K by MDL A simplied

description of the methods is as follows

Initialize the Gaussian mixture model parameters

i

i

W

i

i

K

Repeat

a Estep Compute the model pxj by maximizing the loglikelihood

and assuming the parameter vector to be correct

b Mstep Reevaluate based upon the new data distribution com

puted in the previous step

Until logLj

NK

logK increases by less than

Obtain the Klevel image by encoding in each pixel the cluster mem

bership with the highest likelihood

a

Figure b

c

Figure d

e

Figure f

FIGURE

a Original mammographic image mdb from the MiniMIAS

database b Breast contour and pectoral muscle edge detected

automatically c Eective region of the mammogram obtained after

performing the segmentation steps d Histogram of the eective area of the

mammogram and the mixture of Gaussian components e Fourlevel image

resulting from the EM algorithm f Fibroglandular disc obtained after

thresholding Reproduced with permission from R J Ferrari R M Ran

gayyan R A Borges and A F Fr ere Segmentation of the broglandular

disc in mammograms using Gaussian mixture modelling Medical and

Biological Engineering and Computing

c

IFMBE

Delimit the broglandular disc based upon the density of the pectoral

muscle

Evaluation of the results of segmentation In the work of Ferrari

et al images randomly chosen from the MiniMIAS database

were used to test the segmentation of the broglandular disc All images were

MLO views with m sampling interval and b graylevel quantization

In order to reduce the processing time all images were downsampled with a

xed sampling distance such that the original images with a matrix size of

pixels were reduced to pixels The results obtained

with the downsampled images were mapped to the original mammograms for

subsequent analysis and display

The results obtained were evaluated in consultation with a radiologist ex

perienced in mammography The test images were displayed on a computer

monitor By using the Gimp program the contrast brightness and

zoom options were provided for improved visualization and assessment of the

results of the segmentation procedure

Because the delineation of the broglandular disc is a dicult and time

consuming task and also because the segmentation method may provide dis

joint regions the results were evaluated by using the following subjective pro

cedure The segmented and original images were simultaneously presented to

the radiologist on a computer monitor The radiologist visually compared the

two images and assigned one of ve categories to the result of segmentation

as follows

Excellent Agreement between the segmented disc and the observed disc

on the mammogram is higher than

Good Agreement between the segmented disc and the observed disc on

the mammogram is between and

a

Figure b

c

Figure d

FIGURE

a Breast contour and pectoral muscle edge superimposed on the original

image mdb b Histogram of the eective area of the mammogram

and the mixture of Gaussian components c Threelevel image resulting

from the EM algorithm d Fibroglandular disc obtained after thresholding

Reproduced with permission from R J Ferrari R M Rangayyan R A Borges

and A F Fr ere Segmentation of the broglandular disc in mammograms

using Gaussian mixture modelling Medical and Biological Engineering and

Computing

c

IFMBE

Average Agreement between the segmented disc and the observed disc

on the mammogram is between and

Poor Agreement between the segmented disc and the observed disc on

the mammogram is between and

Complete failure Agreement between the segmented disc and the ob

served disc on the mammogram is less than

The overall results of segmentation of the broglandular disc were consid

ered to be promising Approximately of the cases images were rated

as acceptable for CAD purposes Categories and In Category the

results for images were considered by the radiologist to be underestimated

and the results for images to be overestimated In Category the result

for one image was considered to be underestimated and the results for two

images to be overestimated The results for about of the cases images

in Categories and were considered to be unsatisfactory

In spite of the attractive features of the EM algorithm such as reliable

convergence wellestablished theoretical basis and easy implementation a

major drawback is the problem of convergence to local maxima which is

related to the initialization of the parameters of the model A more ecient

minimization technique such as the deterministic annealing EM algorithm

DAEM may give better performance The DAEM algorithm uses

the principle of maximum entropy and analogy with statistical mechanics to

reformulate the maximization step of the EM algorithm The loglikelihood

function is replaced by a posterior PDF that is parameterized by with

corresponding to the temperature in analogy to the annealing process

This formulation of the PDF provides robustness to DAEM and minimizes

the possibility of convergence to local maxima

The application of the DAEM algorithm to the images in the work of Ferrari

et al provided limited improvement in terms of the t between the

sums of the Gaussians in the mixture models and the corresponding true

histograms However the change in terms of the nal result of segmentation

of the broglandular disc with respect to the result of the EM algorithm was

not signicant

The stochastic EM algorithm is a variant of EM with the aim to avoid

the dependence of the nal results upon the initialization In the stochastic

EM algorithm after the Estep a stochastic classication step the Sstep is

added to obtain partitions of the data according to posterior distributions In

the absence of any other information random initialization is applied The

Mstep uses the partitions to compute new estimates of the mean vector and

the covariance matrix The stochastic EM algorithm presents the following

improvements in comparison with the standard EM algorithm it is adequate

to know an upper bound on the number of classes the solution is essentially

independent of the initialization and the speed of convergence is appreciably

improved

Other features such as texture may also be used as additional information

in order to improve the results Because the segmentation method is based

upon the classication of each individual pixel in the image the resulting

images are often not compact and a followup procedure may be necessary in

order to demarcate the convex hull of the broglandular disc

Motivation for directional analysis of mammograms

Most of the concepts used in image processing and computer vision for ori

ented pattern analysis have their roots in neurophysiological studies of the

mammalian visual system Campbell and Robson suggested that the hu

man visual system decomposes retinal images into a number of ltered images

each of which contains intensity variations over a narrow range of frequency

and orientation Marcelja and Jones and Palmer demonstrated

that simple cells in the primary visual cortex have receptive elds that are

restricted to small regions of space and are highly structured and that their

behavior corresponds to local measurements of frequency

According to Daugman a suitable model for the D receptive

eld proles measured experimentally in mammalian cortical simple cells is

the parameterized family of D Gabor functions see Section Another

important characteristic of Gabor functions or lters is their optimal joint

resolution in both space and frequency which suggests that Gabor lters are

appropriate operators for tasks requiring simultaneous measurement in the

two domains Except for the optimal joint resolution possessed by the

Gabor functions the DoG and dierence of oset Gaussian lters used by

Malik and Perona have similar properties

Gabor lters have been presented in several works on image processing

however most of these works are related to the segmentation and

analysis of texture Rolston and Rangayyan proposed methods for

directional decomposition and analysis of linear components in images using

multiresolution Gabor lters see Section

Multiresolution analysis by using Gabor lters has natural and desirable

properties for the analysis of directional information in images most of these

properties are based upon biological vision studies as described above Other

multiresolution techniques have also been used with success in addressing re

lated topics such as texture analysis and segmentation and image enhance

ment Chang and Kuo for instance developed a new method for texture

classication that uses a treestructured wavelet transform for decomposing

an image Image decomposition was performed by taking into account the en

ergy of each subimage instead of decomposing subsignals in the lowfrequency

channels

Laine and Fan presented a new method for computing features for

texture segmentation based upon a multiscale representation They used a

discrete wavelet packet frame to decompose an image at dierent levels of

resolution The levels were analyzed by using two algorithms for envelope

detection based upon the Hilbert transform and zero crossings Laine and

Fan stated that the selection of the lters for feature extraction should take

into account three important features symmetry frequency response and

boundary accuracy However it should be noted that the Gabor function is

the only function that can achieve optimal joint localization tradeo between

boundary accuracy and good frequency response

Li et al performed analysis and comparison of the directional se

lectivity of wavelet transforms by using three dierent types of frequency

decomposition rectangular hexagonal and radialangular decomposition

The methods were applied to detect spiculated breast masses and lung nod

ules According to Li et al the best results were achieved by using radial

angular decomposition with both mammographic and chest radiographic im

ages Qian et al presented three image processing modules for mass

detection in digital mammography a treestructured lter module for noise

suppression a directional wavelet transform DWT technique for the decom

position of mammographic directional features and a treestructured wavelet

transform for image enhancement By making the parameters of the three

methods adaptive Qian et al improved the results obtained in previous re

lated works In the DWT method which is related to the

technique proposed by Ferrari et al the number of lter orientations

was adaptively selected in the range corresponding to angular band

width in the range

o

o

The adaptive DWT module provided the best

results with an overall classication eciency of

Chang and Laine proposed a method designed with the goal of en

hancing mammograms based upon information regarding orientation They

used a set of separable steerable lters to capture subtle features at dierent

scales and orientations spanning an overcomplete multiscale representation

computed via a fast wavelet transform By using the ltered output images

they generated a coherence image and phase information that were used by

a nonlinear operator applied to the wavelet coecients to enhance the direc

tional information in digital mammograms Finally the enhanced image was

obtained from the modied wavelet coecients via an inverse fast wavelet

transform See Section for further discussion on related topics

Inspired by studies on the Gabor function and the related topics mentioned

above Ferrari et al proposed a scheme based upon a bank of selfsimilar

Gabor functions and the KLT to analyze directional components of images

The method was applied to detect global signs of asymmetry in the broglan

dular discs of the left and right mammograms of a given subject The related

methods are described in the following paragraphs

Directional analysis of broglandular tissue

Ferrari et al used the formulation of D Gabor functions as a Gaussian

modulated by a complex sinusoid 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

Despite this general form there is no standard and precise denition of a D

Gabor function with several variations appearing in the literature

see Sections and Most of these variations are related to the

use of dierent measures of width for the Gaussian envelope and the frequency

of the sinusoid Based upon neurophysiological studies and

wavelet theory the Gabor function normalized in an appropriate

way can be used as a mother wavelet to generate a family of nonorthogonal

Gabor wavelets However as pointed out by Jain and Farrokhnia

although the Gabor function can be an admissible wavelet by removing the

DC response of the function it does not result in an orthogonal decomposi

tion which means that a wavelet transform based upon the Gabor wavelet

includes redundant information A formal mathematical derivation of D Ga

bor wavelets along with the computation of the frame bounds for which this

family of wavelets forms a tight frame is provided by Lee Despite the

lack of orthogonality presented by the Gabor wavelets the Gabor function is

the only function that can achieve the theoretical limit for joint resolution of

information in both the space and the frequency domains

Ferrari et al used the phrase Gabor wavelet representation to rep

resent a bank of Gabor lters normalized to have DC responses equal to zero

and designed in order to have low redundancy in the representation The Ga

bor wavelet representation used was as proposed by Manjunath and Ma

The Gabor wavelets were obtained by dilation and rotation of x y as in

Equation by using the generating function

mn

x y a

m

x

y

a m n integers

x

a

m

x x

cos y y

sin

y

a

m

x x

sin y y

cos

where x

y

is the center of the lter in the spatial domain

n

K

K is

the total number of orientations desired and m and n indicate the scale and

orientation respectively The scale factor a

m

in Equation is meant to

ensure that the energy is independent of m Examples of the Gabor wavelets

used in the work of Ferrari et al are shown in Figure

Equation can be written in the frequency domain as

!u v

u

v

exp

uW

u

v

v

where

u

x

and

v

y

The design strategy used is 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 Fig

ure In this manner it can ensured that the lters will capture most of

the information with minimal redundancy

Other related methods proposed in the literature use either complexvalued

Gabor lters or pairs of Gabor lters with quadraturephase relation

ship In the formulation of Ferrari et al the Gabor wavelet repre

sentation uses only realvalued evensymmetric lters oriented over a range

of

o

only as opposed to the full

o

range commonly described in the

literature Because the Gabor lters are used to extract meaningful features

from real images and hence with Hermitian frequency response the

response to the evensymmetric lter components will remain unchanged for

lters oriented

o

out of phase and the oddsymmetric component will be

negated Thus based upon this fact and also on psychophysical grounds

provided by Malik and Perona Manjunath and Ma ignored one

half of the orientations in their Gabor representation as illustrated in Fig

ure In order to ensure that the bank of Gabor lters designed as above

becomes a family of admissible D Gabor 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 re

sponse at DC This condition was achieved by setting the DC gain of each

lter as ! which ensures that the lters do not respond to regions

with constant intensity

The approach described above results in the following formulas for comput

ing the lter parameters

u

and

v

a

U

h

U

l

S

a b

c d

FIGURE

Examples of Gabor wavelets in the space domain with four orientations

o

o

o

and

o

and four scales

x

and

y

pixels The size of each wavelet image shown is pixels Reproduced

with permission from R J Ferrari R M Rangayyan J E L Desautels and

A F Fr ere Analysis of asymmetry in mammograms via directional ltering

with Gabor wavelets IEEE Transactions on Medical Imaging

c

IEEE

Figure a

b

FIGURE

Examples of Gabor lters in the frequency domain Each ellipse represents

the range of the corresponding lter response from to in squared mag

nitude The plots a and b illustrate two ways of dividing the frequency

spectrum by changing the U

l

U

h

S and K parameters of the Gabor rep

resentation Plot a represents the lter bank used in the work of Ferrari

et al for the analysis of mammograms The redundancy in the rep

resentation is minimized by ensuring that the halfpeak magnitude supports

of the lter responses touch one another Reproduced with permission from

R J Ferrari R M Rangayyan J E L Desautels and A F Fr ere Analysis

of asymmetry in mammograms via directional ltering with Gabor wavelets

IEEE Transactions on Medical Imaging

c

IEEE

u

a U

h

a

p

ln

v

tan

K

h

U

h

u

U

h

ln

i

q

ln

ln

u

U

h

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 frequency of the sinusoid W is set equal to U

h

and m S

Because of the lack of orthogonality of the Gabor wavelets the computation

of the expansion coecients becomes dicult This task however is trivial

when using a set of orthogonal functions because the expansion coecients

given by

c

mn

hfx y

mn

x yi

Z

x

Z

y

fx y

mn

x y dx dy

are the projections of the image fx y onto the same set of functions where

h i denotes the inner product In this case the analysis and synthesis win

dows are the same and the original image can be reconstructed as

fx y

X

m

X

n

hfx y

mn

x yi

mn

x y

However the joint localization and orthogonality of the set of functions are

properties that cannot be met simultaneously Much work has been done to

overcome the problem of the lack of orthogonality of Gabor wavelets with

most of them using dual frames or biorthogonal functions iterative

methods or adjustment of the phasespace sampling in order to ob

tain a reasonably tight frame In the dualframe approach for instance

the set of projection coecients c

mn

D

fx y

e

mn

x y

E

of the dual frame

can be obtained by minimizing the cost function

fx y

X

m

X

n

c

mn

mn

x y

where

e

mn

is the dual frame

In directional ltering and analysis the interest lies in image analysis with

out the requirement of exact reconstruction synthesis of the image There

fore instead of using the wavelet coecients Ferrari et al used the

magnitude of the lter response computed as

a

mn

fx y

even

mn

x y

where

even

mn

x y indicates only the evensymmetric part of the complex Ga

bor lter and represents D convolution

By adjusting the parameters U

l

and U

h

in the Gabor representation of

Manjunath and Ma the range of the frequency spectrum to be used for

multiresolution analysis may be selected The area of each ellipse indicated in

Figure represents the spectrum of frequencies covered by the correspond

ing Gabor lter Once the range of the frequency spectrum is adjusted the

choice of the number of scales and orientation may be made in order to cover

the range of the spectrum as required The choice of the number of scales S

and orientations K used in the work of Ferrari et al for processing

mammographic images was based upon the resolution required for detecting

oriented information with high selectivity The frequency band

widths of the simple and complex cells in the mammalian visual system 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 pro

cessing mammographic images Ferrari et al indirectly set the Gabor

representation to have a frequency bandwidth of approximately one octave

and an angular bandwidth of

o

The eects of changing the U

l

U

h

S and

K parameters of the Gabor representation as above on frequency localization

are shown in Figure

The directional analysis method proposed by Ferrari et al starts by

computing the Gabor wavelets using four scales S and twelve directions

K with U

l

and U

h

cycles pixel The parameters U

l

and

U

h

were chosen according to the scales of the details of interest in the mam

mographic images Diering from other Gabor representations

it should be noted that the parameters U

l

U

h

S and K in the representation

described above have to be adjusted in a coordinated manner by taking into

account the desirable frequency and orientation bandwidths In the Gabor

representation of Manjunath and Ma the lters are designed so as to

represent an image with minimal redundancy the ellipses as in Figure

touch one another

The lter outputs for each orientation and the four scales were analyzed

by using the KLT see Section The KLT was used to select the

principal components of the lter outputs preserving only the most relevant

directional elements present at all of the scales considered Results were then

combined as illustrated in Figure in order to allow the formation of an

Sdimensional vector x for each pixel from each set of the corresponding

pixels in the ltered images S

The vectors corresponding to each position in the lter responses were used

to compute the mean vector and the covariance matrix The eigenvalues

and eigenvectors of the covariance matrix were then computed and arranged in

descending order in a matrixA such that the rst row ofA was the eigenvector

corresponding to the largest eigenvalue and the last row was the eigenvector

corresponding to the smallest eigenvalue The rst N principal components

corresponding to of the total variance were then selected and used to

represent the oriented components at each specic orientation The principal

FIGURE

Formation of the vector x x from the corresponding pixels of the same

orientation and four scales Reproduced with permission from R J Ferrari

R M Rangayyan J E L Desautels and A F Fr ere Analysis of asymmetry

in mammograms via directional ltering with Gabor wavelets IEEE Trans

actions on Medical Imaging

c

IEEE

components were computed as y A x Analysis of variance was

performed by evaluating the eigenvalues of the matrix A The KLT method

is optimal in the sense that it minimizes the MSE between the vectors x and

their resulting approximations y The result of application of the KLT to all

orientations as described above is a set of K images where K is the number

of orientations

Because Gabor wavelets are nonorthogonal functions they do not have a

perfect reconstruction condition This fact results in a small amount of out

ofband energy interfering with the reconstruction which is translated into

artifacts in the reconstructed image Although reconstruction of the original

image from the ltered images is not required in directional analysis the

eects of spectral leakage need to be removed For this purpose the images

resulting from the KLT were thresholded by using the maximum of Otsus

threshold value see Section computed for the K images

Phase and magnitude images indicating the local orientation and inten

sity were composed by vector summation of the K ltered images as

illustrated in Figure Rose diagrams were composed from the phase and

magnitude images to represent the directional distribution of the broglan

dular tissue in each mammogram The complete algorithm for directional

analysis based upon Gabor ltering is summarized in Figure

FIGURE

Block diagram of the procedure for directional analysis using Gabor wavelets

Reproduced with permission from R J Ferrari R M Rangayyan J E L De

sautels and A F Fr ere Analysis of asymmetry in mammograms via direc

tional ltering with Gabor wavelets IEEE Transactions on Medical Imaging

c

IEEE

Characterization of bilateral asymmetry

Figures and show two pairs of images from the MiniMIAS database

All of the images in this database are pixels in size with

m sampling interval and b graylevel quantization Figures a

and b are respectively the MIAS images mdb and mdb representing

the right and left MLO mammograms of a woman classied as a normal

case Figures a and b show the MIAS images mdb and mdb

classied as a case of architectural distortion

Only the broglandular disc of each mammogram was used to compute the

directional components due to the fact that most of the directional compo

nents such as connective tissues and ligaments exist in this specic region of

the breast The broglandular disc ROIs of the mammograms selected for

illustration are shown in Figures c d c and d

Figure shows the principal components obtained by applying the KLT

to the Gabor lter outputs at the orientation of

o

and four scales for

the ROI in Figure d It can be seen that the relevant information is

concentrated in the rst two principal components This is evident based

upon an evaluation of the eigenvalues listed in the caption of Figure In

this example only the rst two principal components were used to represent

the oriented components in the

o

orientation because their eigenvalues

add to of the total variance After thresholding the ltered

images with Otsus method in order to eliminate the eects of spectral leakage

the magnitude and phase images were composed by vector summation as

illustrated in Figures and

Figures a d show the result of Gabor ltering for the normal

case in Figure The rose diagrams in Figures c and d show the

distribution of the tissues in the broglandular discs of both the left and right

views An inspection of the rose diagrams shows that the results obtained

are in good agreement with visual analysis of the ltered results in Figures

a and b and the corresponding ROIs in Figures c and d The

most relevant angular information indicated in the rose diagrams are similar

The results of the ltering process for the case of architectural distortion

see Figure along with the respective rose diagrams are shown in Fig

ure By analyzing the results of ltering we can notice a modication of

the tissue pattern caused by the presence of a highdensity region An impor

tant characteristic of the Gabor lters may be seen in the result the lters

do not respond to regions with nearly uniform intensity that is to regions

without directional information see Figures c and a This is

an important outcome that could be used to detect asymmetric dense regions

and local foci of architecture distortion The global distortion of the tissue

ow pattern is readily seen by comparison of the rose diagrams of the left and

right breasts

a b

c d

FIGURE

Images mdb and mdb of a normal case a and b Original im

ages pixels at m pixel The breast boundary white and

pectoral muscle edge black detected are also shown c and d Fibroglan

dular discs segmented and enlarged pixels Histogram equalization

was applied to enhance the global contrast of each ROI for display purposes

only Reproduced with permission from R J Ferrari R M Rangayyan J E L

Desautels and A F Fr ere Analysis of asymmetry in mammograms via direc

tional ltering with Gabor wavelets IEEE Transactions on Medical Imaging

c

IEEE

a b

c d

FIGURE

Images mdb and mdb of a case of architectural distortion a

and b Original images pixels at m pixel The breast

boundary white and pectoral muscle edge black detected are also shown

c and d Fibroglandular discs segmented and enlarged pixels

Histogram equalization was applied to enhance the global contrast of each

ROI for display purposes only Reproduced with permission from R J Fer

rari R M Rangayyan J E L Desautels and A F Fr ere Analysis of asym

metry in mammograms via directional ltering with Gabor wavelets IEEE

Transactions on Medical Imaging

c

IEEE

a b

c d

FIGURE

The images a b c and d are respectively the rst second third

and fourth components resulting from the KLT applied to the Gabor lter re

sponses with orientation

o

to the ROI of the image mdb shown in Fig

ure d The eigenvalues of the four components above are

and

The images were full brightnesscontrast

corrected for improved visualization Reproduced with permission from R J

Ferrari R M Rangayyan J E L Desautels and A F Fr ere Analysis of

asymmetry in mammograms via directional ltering with Gabor wavelets

IEEE Transactions on Medical Imaging

c

IEEE

The rose diagrams in Figures and present a strong visual associ

ation with the directional components of the corresponding images and may

be used by radiologists as an aid in the interpretation of mammograms

Feature extraction and pattern classi cation In order to characterize

bilateral asymmetry in an objective manner three features were derived the

entropy H Equation the rst moment M

Equation and the

second central moment or variance M

Equation of the rose diagram

given by the dierence between the rose diagrams computed for the left and

right mammograms of the same individual

Classication of the normal and asymmetric cases was conducted by using

the Bayesian linear classier see Section The Gaussian dis

tribution was assumed in order to model the PDF and the parameters of

the model were estimated by using the training samples The prior proba

bilities of the normal and asymmetry classes were assumed to be equal and

the covariance matrix was calculated in a pooled manner by averaging the

covariance matrices of the normal and asymmetric classes The leaveoneout

methodology was used to estimate the classication accuracy

The directional analysis scheme was applied to images normal cases

cases of asymmetry and six cases of architectural distortion from the

MiniMIAS database An exhaustive combination approach was used to

select the best set of features The selection was conducted based upon the

classication results obtained by using the leaveoneout method The best

result by using only one feature in the classication process was achieved

by the rstorder angular moment M

with the sensitivity specicity and

average accuracy values equal to and respectively When

using two features the best result was achieved with the combination of the

rstorder angular moment M

and the entropy H features indicating

that of the asymmetric and distortion cases and of the normal

cases were correctly classied The average rate of correct classication in

this case was The low rate of specicity may be explained by the

fact that even normal cases present natural signs of mild asymmetry the

mammographic imaging procedure may also distort the left and right breasts

in dierent ways

In a subsequent study Rangayyan et al revised the directional analysis

procedures as shown in Figure The rose diagrams of the left and right

mammograms were aligned such that their mean angles corresponded to the

straight line perpendicular to the pectoral muscle and then subtracted to

obtain the dierence rose diagram In addition to the features H M

and

M

of the dierence rose diagram as described above the dominant orientation

R

and circular variance s

were computed as follows

X

R

N

X

i

R

i

cos

i

a b

c d

FIGURE

Results obtained for the normal case in Figure a and b Magnitude

images c and d Rose diagrams The magnitude images were histogram

equalized for improved visualization The rose diagrams have been congured

to match the mammograms in orientation Reproduced with permission from

R J Ferrari R M Rangayyan J E L Desautels and A F Fr ere Analysis

of asymmetry in mammograms via directional ltering with Gabor wavelets

IEEE Transactions on Medical Imaging

c

IEEE

a b

c d

FIGURE

Results obtained for the case of architectural distortion in Figure a and

b Magnitude images c and d Rose diagrams The magnitude images

were histogramequalized for improved visualization The rose diagrams have

been congured to match the mammograms in orientation Reproduced with

permission from R J Ferrari R M Rangayyan J E L Desautels and A F

Fr ere Analysis of asymmetry in mammograms via directional ltering with

Gabor wavelets IEEE Transactions on Medical Imaging

c

IEEE

Y

R

N

X

i

R

i

sin

i

R

atan

Y

R

X

R

and

s

q

X

R

Y

R

where R

i

is the normalized value and

i

is the central angle of the i

th

angle

band of the dierence rose diagram and N is the number of bins in the rose

diagram

In addition a set of features including seven of Hus moments see Sec

tion and Equation and the area average density eccentricity and

stretch were computed to characterize the shape of the segmented broglan

dular discs Eccentricity was computed as

m

m

m

m

m

where m

pq

are the geometric invariant moments as described in Section

The stretch parameter was computed as

x

max

x

min

y

max

y

min

where x

max

x

min

y

max

and y

min

are the corner coordinates of the rectangle

delimiting the broglandular disc Feature selection was performed by PCA

and exhaustive combination techniques With PCA only the components

associated with of the total variance were used in the classication step

Classication was performed using linear and quadratic Bayesian classiers

with the leaveoneout method

The revised directional analysis scheme was applied to images nor

mal cases cases of asymmetry and eight cases of architectural distortion

from the MiniMIAS database The best overall classication accuracy

of with a sensitivity of and specicity of was obtained us

ing the four features

R

M

M

and H computed from the aligneddierence

rose diagrams using the quadratic classier The morphometric measures and

moments after PCAbased feature selection resulted in an overall classica

tion accuracy of only with the linear classier The combination of

all of the directional statistics morphometric measures and moments after

PCAbased feature selection resulted in an overall classication accuracy of

with a sensitivity of and specicity of with the linear

classier The results indicate the importance of directional analysis of the

broglandular tissue in the detection of bilateral asymmetry

FIGURE

Block diagram of the revised procedure for the analysis of bilateral asymme

try

Application Architectural Distortion in Mammo

grams

Mammography is the best available tool for detecting early breast cancer

screening programs have been shown to reduce mortality rates by to

Chapter However the sensitivity of screening mammography

is aected by image quality and the radiologists level of expertise Bird et

al estimated the sensitivity of screening mammography to be between

and they also observed that misinterpretation of breast cancer signs

accounts for of the errors and overlooking signs is responsible for of

the missed lesions The extent of errors due to overlooking of lesions reinforces

the need for CAD tools in mammography

CAD techniques and systems have been proposed to enhance the sensitiv

ity of the detection of breast cancer although these techniques are eective

in detecting masses and calcications they have been found to fail in the

detection of architectural distortion with adequate levels of accuracy

Therefore new CAD systems are desirable for targeted detection of subtle

mammographic abnormalities such as architectural distortion which are the

most frequent source of screening errors and falsenegative ndings related to

cases of interval cancer

Architectural distortion is dened in BIRADS

TM

as follows The

normal architecture of the breast is distorted with no denite mass visible

This includes spiculations radiating from a point and focal retraction or dis

tortion at the edge of the parenchyma Focal retraction is considered to be

easier to perceive than spiculated distortion within the breast parenchyma

p Architectural distortion could be categorized as malignant or benign

the former including cancer and the latter including scar and softtissue dam

age due to trauma

According to van Dijck et al in nearly half of the screendetected

cancers minimal signs appeared to be present on the previous screening mam

mogram two years before the diagnosis Burrell et al in a study of

screening interval breast cancers showed that architectural distortion is the

most commonly missed abnormality in falsenegative cases Sickles re

ported that indirect signs of malignancy such as architectural distortion

bilateral asymmetry single dilated duct and developing densities account

for almost of the detected cancers Broeders et al suggested that

improvement in the detection of architectural distortion could lead to an ef

fective improvement in the prognosis of breast cancer patients

Detection of spiculated lesions and distortion

The breast contains several piecewise linear structures such as ligaments

ducts and blood vessels that cause oriented texture in mammograms The

presence of architectural distortion is expected to change the normal oriented

texture of the breast

Figure a shows a mammogram of a normal breast where the normal

oriented texture is observed The mammogram in Figure b exhibits

architectural distortion Figure shows an enlarged view of the site of

architectural distortion Observe the appearance of lines radiating from a

central point in the ROI It is also noticeable that the perceived increased

density close to the center of the ROI is due to the broglandular disk as it

can be observed in the full mammogram view in Figure b

a b

FIGURE

a Mammogram showing a normal breast image mdb from the Mini

MIAS database Width of image pixels mm b Archi

tectural distortion present in a mammogram from the MiniMIAS database

mdb Width of image pixels mm The square box overlaid

on the gure represents the ROI including the site of architectural distortion

shown enlarged in Figure Reproduced with permission from F J Ayres

and R M Rangayyan Characterization of architectural distortion in mam

mograms via analysis of oriented texture IEEE Engineering in Medicine and

Biology Magazine January

c

IEEE

FIGURE

Detail of mammogram mdb showing the site of architectural distortion

marked by the box in Figure b Width of image pixels mm

Reproduced with permission from F J Ayres and R M Rangayyan Charac

terization of architectural distortion in mammograms via analysis of oriented

texture IEEE Engineering in Medicine and Biology Magazine January

c

IEEE

Several researchers have directed attention to the detection of stellate and

spiculated patterns associated with masses in mammograms see Mudigonda

et al and Section for a review on this subject However most of such

methods have been directed toward and applied to the detection of tumors

with spiculations A common strategy in most of the methods reported in

the literature for this application is to assume a stellate appearance for the

lesion detection is performed by tting a star pattern to the texture of the

breast parenchyma

Qian et al proposed a directional wavelet transform to detect spicules

radiating from masses Kegelmeyer et al used local edge orientation his

tograms and Laws texture features to detect stellate lesions A sensitivity of

was obtained with a specicity of Karssemeijer and te Brake

applied operators sensitive to radial patterns of straight lines to an orienta

tion eld to detect stellate patterns and obtained a sensitivity of with

one false positive per image however only lesions with radiating spicules re

gardless of the presence of a central density are detected in their method

Mudigonda et al investigated the detection of masses using density slic

ing and texture oweld a sensitivity of with false positives per

image was achieved see Section

Zwiggelaar et al proposed a technique to detect abnormal patterns of

linear structures by detecting the radiating spicules and"or the central mass

expected to occur with spiculated lesions PCA was applied to a training

set of mammograms including normal tissue patterns and spiculated lesions

The results of PCA were used to construct a basis set of oriented texture

patterns which was used to analyze radiating structures A sensitivity of

was obtained in the classication of pixels as belonging to normal tissue

or lesions but with low specicity Sampat and Bovik employed ltering

in the Radon domain to enhance mammograms followed by the usage of radial

spiculation lters to detect spiculated lesions however statistical evaluation

of the performance of the technique was not conducted

Mudigonda and Rangayyan proposed the use of texture oweld to

detect architectural distortion based upon the local coherence of texture ori

entation only preliminary results were given indicating the potential of the

technique in the detection of architectural distortion Ferrari et al used

Gabor lters to detect oriented patterns with subsequent analysis designed

to characterize global bilateral asymmetry however the method was not

aimed at detecting focal architectural distortion see Section Matsubara

et al used mathematical morphology to detect architectural distortion

around the skin line and a concentration index to detect architectural dis

tortion within the mammary gland they reported sensitivity rates of

with false positives per image and with false positives per image

respectively in the two tasks

Burhenne et al studied the performance of a commercial CAD sys

tem in the detection of masses and calcications in screening mammography

obtaining a sensitivity of in the detection of masses and architectural

distortion Evans et al investigated the ability of a commercial CAD

system to mark invasive lobular carcinoma of the breast the system identied

correctly of cases of architectural distortion Birdwell et al eval

uated the performance of a commercial CAD system in marking cancers that

were overlooked by radiologists the software detected ve out of six cases

of architectural distortion However Baker et al found the sensitivity

of two commercial CAD systems to be poor in detecting architectural distor

tion fewer than of the cases of architectural distortion were detected

These ndings indicate the need for further research in this area and the

development of algorithms designed specically to characterize architectural

distortion

Ayres and Rangayyan proposed the application of Gabor

lters and phase portraits to characterize architectural distortion in ROIs se

lected from mammograms as well as to detect sites of architectural distortion

in full mammograms their methods and results are described in the following

subsections

Phase portraits

Phase portraits provide an analytical tool to study systems of rstorder dif

ferential equations The method has proved to be useful in characterizing

oriented texture

Let pt and qt denote two dierentiable functions of time t related by a

system of rstorder dierential equations as

#pt F pt qt

#qt Gpt qt

where the dot above the variable indicates the rstorder derivative of the

function with respect to time and F and G represent functions of p and q

Given initial conditions p and q the solution pt qt to Equation

can be viewed as a parametric trajectory of a hypothetical particle in the

p q plane placed at p q at time t and moving through the p q

plane with velocity #pt #qt The p q plane is referred to as the phase

plane of the system of rstorder dierential equations The path traced by

the hypothetical particle is called a streamline of the vector eld #p #q The

phase portrait is a graph of the possible streamlines in the phase plane A

xed point of Equation is a point in the phase plane where #pt and

#qt a particle left at a xed point remains stationary

When the system of rstorder dierential equations is linear Equation

assumes the form

#pt

#qt

A

pt

qt

b

where A is a matrix and b is a column matrix a vector In

this case there are only three types of phase portraits node saddle and

spiral The type of phase portrait can be determined from the nature of

the eigenvalues of A as shown in Table The center p

q

of the phase

portrait is given by the xed point of Equation

#pt

#qt

p

q

A

b

Solving Equation yields a linear combination of complex exponentials

for pt and qt whose exponents are given by the eigenvalues of A multiplied

by the time variable t Table illustrates the streamlines obtained by solving

Equation for a node a saddle and a spiral phase portrait the solid lines

indicate the movement of the pt and the qt components of the solution

and the dashed lines indicate the streamlines The formation of each phase

portrait type is explained as follows

Node the components pt and qt are exponentials that either simul

taneously converge to or diverge from the xedpoint coordinates p

and q

Saddle the components pt and qt are exponentials while one of the

components either pt or qt converges to the xed point the other

diverges from the xed point

Spiral the components pt and qt are exponentially modulated sinu

soidal functions the resulting streamline forms a spiral curve

Associating the functions pt and qt with the x and y coordinates of

the Cartesian image plane we can dene the orientation eld generated by

Equation as

x yjAb arctan

#qt

#pt

which is the angle of the velocity vector #pt #qt with the x axis at x y

pt qt Table lists the three phase portraits and the corresponding

orientation elds generated by a system of linear rstorder dierential equa

tions

Using the concepts presented above the orientation eld of a textured image

may be described qualitatively by determining the type of the phase portrait

that is most similar to the orientation eld along with the center of the phase

portrait This notion was employed by Ayres and Rangayyan

to characterize architectural distortion

Estimating the orientation eld

The analysis of oriented texture is an important task in computer vision

and several methods have been proposed for this task such as weighted an

gle histograms integral curves and phase portraits

Ayres and Rangayyan analyzed the orientation eld in ROIs of

mammograms through the usage of phase portraits in order to characterize

architectural distortion In order to extract the texture orientation at each

pixel of the ROI the ROI was ltered with a bank of Gabor lters of dierent

orientations see Section The Gabor lter kernel oriented at the

angle was formulated as

gx y

x

y

exp

x

x

y

y

cosf

o

x

Kernels at other angles were obtained by rotating this kernel A set of

kernels was used with angles spaced evenly over the range

Gabor lters may be used as line detectors In the work of Ayres and

Rangayyan the parameters in Equation namely

x

y

and f

o

were derived from a design rule as follows Let be the thickness of

the line detector This parameter constrains

x

and f

o

as follows

The amplitude of the exponential term in Equation that is the

Gaussian term is reduced to one half of its maximum at x and

y therefore

x

p

ln

TABLE

Phase Portraits for a System of Linear Firstorder Dierential

Equations

Phase

portrait

type

Eigenvalues

Streamlines

Appearance of the

orientation eld

Node

Real eigenvalues

of same sign

Saddle

Real eigenvalues

of opposite sign

Spiral

Complex

eigenvalues

Solid lines indicate the movement of the pt and the qt components of the

solution dashed lines indicate the streamlines Reproduced with permission

from F J Ayres and R M Rangayyan Characterization of architectural dis

tortion in mammograms via analysis of oriented texture IEEE Engineering

in Medicine and Biology Magazine January

c

IEEE

The cosine term has a period of therefore f

o

The value of

y

was dened as

y

l

x

where l determines the elongation

of the Gabor lter in the orientation direction with respect to its thickness

The values pixels corresponding to a thickness of mm at a pixel size

of m and l were determined empirically by observing the typical

spicule width and length in mammograms with architectural distortion in the

MiniMIAS database

The eects of the dierent design parameters are shown in Figure and

are as follows

Figures a and e show the impulse response of a Gabor lter and

its Fourier transform magnitude respectively

In Figure b the Gabor lter of Figure a is stretched in

the x direction by increasing the elongation factor l Observe that the

Fourier spectrum of the new Gabor lter shown in Figure f is

compressed in the horizontal direction

The Gabor lter shown in Figure c was obtained by increasing

the parameter of the original Gabor lter thus enlarging the lter in

both the x and y directions Correspondingly the Fourier spectrum of

the enlarged lter shown in Figure g has been shrunk in both

the vertical and horizontal directions

The eect of rotating the Gabor lter by

o

counterclockwise is dis

played in Figures d and h that show the rotated Gabor lters

impulse response and the corresponding Fourier spectrum

The texture orientation at a pixel was estimated as the orientation of the

Gabor lter that yielded the highest magnitude response at that pixel The

orientation at every pixel was used to compose the orientation eld The

magnitude of the corresponding lter response was used to form themagnitude

image The magnitude image was not used in the estimation of the phase

portrait but was found to be useful for illustrative purposes

Let x y be the texture orientation at x y and g

k

x y k

be the Gabor lter oriented at

k

k Let fx y be the ROI

of the mammogram being processed and f

k

x y f g

k

x y where the

asterisk denotes linear D convolution Then the orientation eld of fx y

is given by

x y

k

max

where k

max

argfmax

k

jf

k

x yjg

Characterizing orientation elds with phase portraits

In the work of Ayres and Rangayyan the analysis of oriented

texture patterns was performed in a twostep process First the orientation

a b c d

e f g h

FIGURE

Eects of the dierent parameters of the Gabor lter a Example of the

impulse response of a Gabor lter b The parameter l is increased the

Gabor lter is elongated in the x direction c The parameter is increased

the Gabor lter is enlarged in the x and y directions d The angle of the

Gabor lter is modied Figures e h correspond to the magnitude of the

Fourier transforms of the Gabor lters in a d respectively The fre

quency component is at the center of the spectra displayed Reproduced with

permission from F J Ayres and R M Rangayyan Characterization of archi

tectural distortion in mammograms via analysis of oriented texture IEEE

Engineering in Medicine and Biology Magazine January

c

IEEE

eld x y of the ROI was computed in a small analysis window The sliding

analysis window was centered at pixels within the ROI avoiding window posi

tions with incomplete data at the edges of the ROI for the estimation ofA and

b Second the matrix A and the vector b in Equation were estimated

such that the best match was achieved between x y and x yjAb The

eigenvalues of A determine the type of the phase portrait present in x y

the xed point of the phase portrait is given by Equation

The estimates of A and b were obtained as follows For every point x y

let $x y sinx y x yjAb represent the error between the orien

tation of the texture given by Equation and the orientation of the model

given by Equation Then the estimation problem is that of nding A

and b that minimize the sum of the squared error

X

x

X

y

$

x y

X

x

X

y

fsinx y x yjAbg

which may be solved using a nonlinear leastsquares algorithm

The ROI was investigated by sliding the analysis window through the ori

entation eld of the ROI and accumulating the information obtained that

is the type of the phase portrait and the location of the xed point for each

window position as follows

Create three maps one for each type of phase portrait hereafter called

the phase portrait maps that will be used to accumulate information

from the sliding analysis window The maps are initialized to zero and

are of the same size as the ROI or the image being processed

Slide the analysis window through the orientation eld of the ROI At

each position of the sliding window determine the type of the phase

portrait and compute the xed point of the orientation eld

Increment the value at the location of the xed point in the correspond

ing phase portrait map

The size of the sliding analysis window was set at pixels mm

The three maps obtained as above provide the results of a voting procedure

and indicate the possible locations of xed points corresponding to texture

patterns that approximately match the node saddle and spiral phase por

traits It is possible that for some positions of the sliding analysis window

the location of the xed point falls outside the spatial limits of the ROI or

image being processed the votes related to such results were ignored The

value at each location x y in a phase portrait map provides the degree of

condence in determining the existence of the corresponding phase portrait

type centered at x y The three phase portraits were expected to relate to

dierent types of architectural distortion

Feature extraction for pattern classication

The estimates of the xedpoint location for a given phase portrait pattern can

be scattered around the true xedpoint position due to the limited precision

of the estimation procedure the presence of multiple overlapping patterns

the availability of limited data within the sliding analysis window and the

presence of noise A local accumulation of the votes is necessary to diminish

the eect of xedpoint location errors Ayres and Rangayyan

employed a Gaussian smoothing lter with a standard deviation of pixels

mm for this purpose

For the purpose of pattern classication six features were extracted to char

acterize each ROI the maximum of each phase portrait map three features

and the entropy of each phase portrait map three features The maximum

of each map conveys information about the likelihood of the presence of the

corresponding phase portrait type and the entropy relates to the uncertainty

in the location of the xed point in each map The entropy H of a map hx y

was computed as

H hx y

X

x

X

y

hx y

S

h

ln

hx y

S

h

where

S

h

X

x

X

y

hx y

A map with a dense spatial concentration of votes is expected to have a large

maximum value and a low entropy On the contrary a map with a wide

scatter of votes may be expected to have a low maximum and a large entropy

Application to segments of mammograms

Ayres and Rangayyan analyzed a set of ROIs each of size

pixels mm with a resolution of m selected from

the MiniMIAS database The set included ROIs with architec

tural distortion all the cases of architectural distortion available in the MIAS

database ROIs with normal tissue patterns eight ROIs with spiculated

malignant masses four ROIs with circumscribed malignant masses ROIs

with spiculated benign masses ROIs with circumscribed benign masses

and two ROIs with malignant calcications The size of the ROIs was chosen

to accommodate the largest area of architectural distortion or mass identi

ed in the MiniMIAS database ROIs related to all of the masses in the

database were included The normal ROIs included examples of overlapping

ducts ligaments and other parenchymal patterns Only the central portion of

pixels of each ROI was investigated using a sliding analysis window

of size pixels the remaining outer ribbon of pixels was not processed

in order to discard the eects of the preceding ltering steps

TABLE

Results of Linear Discriminant Analysis for ROIs with

Architectural Distortion Using the Leaveoneout Method

Architectural %ROIs Classied as

distortion Architectural distortion Other

Benign

Malignant

Total TP FN

TP true positives FN false negatives The results correspond to the

prior probability of belonging to the architectural distortion class being

Sensitivity Reproduced with permission from F J Ayres and R M

Rangayyan Characterization of architectural distortion in mammograms via

analysis of oriented texture IEEE Engineering in Medicine and Biology Mag

azine January

c

IEEE

Figure illustrates the results obtained for an image with architectural

distortion mdb The maximum of the node map is larger than the max

ima of the other two maps Also the scattering of votes in the node map is

less than that in the saddle and spiral maps These results indicate that the

degree of scattering of the votes quantied by the entropy of the correspond

ing map and the maximum of each of the three phase portrait maps could

be useful features to distinguish between architectural distortion and other

patterns

Linear discriminant analysis was performed using SPSS with stepwise

feature selection Architectural distortion was considered as a positive nding

with all other test patterns normal tissue masses and calcications being

considered as negative ndings The statistically signicant features were

the entropy of the node map and the entropy of the spiral map the other

features were deemed to be not signicant by the statistical analysis package

and were discarded in all subsequent analysis With the prior probability of

architectural distortion set to the sensitivity obtained was and

the specicity was The fraction of cases correctly classied was

Tables and present the classication results with the prior probability

of architectural distortion being An overall classication accuracy of

was achieved

Detection of sites of architectural distortion

Ayres and Rangayyan hypothesized that architectural distor

tion would appear as an oriented texture pattern that can be locally approxi

mated by a linear phase portrait model furthermore it was expected that the

a b

c

d e f

FIGURE

Analysis of the ROI from the image mdb which includes architectural dis

tortion a ROI of size pixels mm b magnitude image

c orientation eld superimposed on the original ROI d node map with

intensities mapped from to e saddle map mapped to

f spiral map mapped to This image was correctly

classied as belonging to the architectural distortion category Table

Reproduced with permission from F J Ayres and R M Rangayyan Charac

terization of architectural distortion in mammograms via analysis of oriented

texture IEEE Engineering in Medicine and Biology Magazine January

c

IEEE

TABLE

Results of Linear Discriminant Analysis for ROIs Without

Architectural Distortion Using the Leaveoneout Method

Type %ROIs Classied as

Architectural distortion Other

CB

Masses SB

CM

SM

Calcications

Normal

Total FP TN

CB circumscribed benign mass CM circumscribed malignant tumor SB

spiculated benign mass SM spiculated malignant tumor FP false pos

itives TN true negatives The results correspond to the prior probability

of belonging to the architectural distortion class being Specicity

Reproduced with permission from F J Ayres and R M Rangayyan

Characterization of architectural distortion in mammograms via analysis of

oriented texture IEEE Engineering in Medicine and Biology Magazine Jan

uary

c

IEEE

xedpoint location of the phase portrait model would fall within the breast

area in the mammogram Then the numbers of votes cast at each position of

the three phase portrait maps would indicate the likelihood that the position

considered is a xed point of a node a saddle or a spiral pattern

Before searching the maps for sites of distortion the orientation eld was

ltered and downsampled as follows Let hx y be a Gaussian lter of stan

dard deviation

h

dened as

hx y

h

exp

x

y

h

Dene the images sx y sinx y and cx y cosx y where

x y is the orientation eld Then the ltered orientation eld

f

x y is

obtained as

f

x y

arctan

hx y sx y

hx y cx y

where the asterisk denotes D convolution

The ltered orientation eld was downsampled by a factor of four thus

producing the downsampled orientation eld

d

as

d

x y

f

x y

The ltering and downsampling procedures summarized in Figure

were applied in order to reduce noise and and also to reduce the computational

eort required for the processing of full mammograms The ltering procedure

described above is a variant of Raos dominant local orientation method

a Gaussian lter has been used instead of a box lter

The following procedure was used to detect and locate sites of architectural

distortion using only the node map

The node map is ltered with a Gaussian lter of standard deviation

equal to pixel mm

The ltered node map is thresholded with the same threshold value for

all images

The thresholded image is subjected to the following series of morpholog

ical operations to group positive responses that are close to one another

and to reduce each region of positive response to a single point The re

sulting points indicate the detected locations of architectural distortion

a A closing operation is performed to group clusters of points that

are less than mm apart The structural element is a disk of radius

pixels mm

b A close holes lter is applied to the image The resulting image

includes only compact regions

FIGURE

Filtering and downsampling of the orientation eld Figure courtesy of F J

Ayres

c The image is subjected to a shrink lter where each compact

region is shrunk to a single pixel

The threshold value inuences the sensitivity of the method and the number

of false positives per image A high threshold value reduces the number of

false positives but also reduces the sensitivity A low threshold value increases

the number of false positives

The method was applied to mammograms exhibiting architectural dis

tortion selected from the MiniMIAS database The mammograms were

MLO views digitized to pixels at a resolution of m and

b pixel Figures and illustrate the steps of the method as applied

to image mdb Observe that the ltered orientation eld Figure d

is smoother and more coherent as compared to the original orientation eld

Figure c the pattern of architectural distortion is displayed better in

the ltered orientation eld

The architectural distortion present in the mammogram mdb has a stel

late or spiculated appearance As a consequence a large number of votes

have been cast into the node map at a location close to the center of the

distortion as seen in Figure c Another point of accumulation of votes

in the node map is observed in Figure c at the location of the nipple

This is not unexpected the breast has a set of ducts that carry milk to the

nipple the ducts appear in mammograms as linear structures converging to

the nipple Observe that the node map has the strongest response of all maps

Figure a b

within the site of architectural distortion given by the MiniMIAS database

The technique has resulted in the identication of two locations as sites of

architectural distortion one true positive and one false positive as shown in

Figure d

The freeresponse receiver operating characteristics FROC was derived

by varying the threshold level in the detection step the result is shown in

Figure See Section for details on ROC analysis A sensitivity

of was obtained at false positives per image Further work is required

in order to reduce the number of false positives and improve the accuracy of

detection

Remarks

Preferred orientation and directional distributions relate to the functional

integrity of several types of tissues and organs changes in such patterns could

indicate structural damage as well as recovery Directional analysis could

c d

FIGURE

a Image mdb from the MiniMIAS database The circle indicates

the location and the extent of architectural distortion as provided in the Mini

MIAS database b Magnitude image after Gabor ltering c Orien

tation eld superimposed on the original image Needles have been drawn for

every fth pixel d Filtered orientation eld superimposed on the original

image Reproduced with permission from F J Ayres and R M Rangayyan

Detection of architectural distortion in mammograms using phase portraits

Proceedings of SPIE Medical Imaging Image Processing Volume

pp

c

SPIE See also Figure

Figure a b

c d

FIGURE

Phase portrait maps derived from the orientation eld in Figure d and

the detection of architectural distortion a Saddle map values are scaled

from the range to b Spiral map values are scaled from the

range to a Node map values are scaled from the range

to d Detected sites of architectural distortion superimposed

on the original image the solid line indicates the location and spatial extent

of architectural distortion as given by the MiniMIAS database the

dashed lines indicate the detected sites of architectural distortion one true

positive and one false positive Reproduced with permission from F J Ayres

and R M Rangayyan Detection of architectural distortion in mammograms

using phase portraits Proceedings of SPIE Medical Imaging Image

Processing Volume pp

c

SPIE

FIGURE

Freeresponse receiver operating characteristics FROC curve for the detec

tion of sites of architectural distortion Reproduced with permission from F J

Ayres and R M Rangayyan Detection of architectural distortion in mam

mograms using phase portraits Proceedings of SPIE Medical Imaging

Image Processing Volume pp

c

SPIE

therefore be used to study the health and wellbeing of a tissue or organ as

well as to follow the pathological and physiological processes related to injury

treatment and healing

In this chapter we explored the directional characteristics of several biomed

ical images We have seen several examples of the application of fan lters

and Gabor wavelets The importance of multiscale or multiresolution analysis

in accounting for variations in the size of the pattern elements of interest has

been demonstrated In spite of theoretical limitations the methods for direc

tional analysis presented in this chapter have been shown to lead to practically

useful results in important applications

Study Questions and Problems

Selected data les related to some of the problems and exercises are available at the

site

wwwenelucalgarycaPeopleRangaenel

Discuss how entropy can characterize a directional distribution

Discuss the limitations of fan lters Describe how Gabor functions address

these limitations

Using an image with line segments of various widths as an example discuss

the need for multiscale or multiresolution analysis

Laboratory Exercises and Projects

Prepare a test image with line segments of di erent directions lengths and

widths with overlap Apply Gabor lters at a few di erent scales and angles

as appropriate Evaluate the results in terms of

a the lengths of the extracted components and

b the widths of the extracted components

Discuss the limitations of the methods and the artifacts in the results

Decompose your test image in the preceding problem using eight sector or fan

lters spanning the full range of

o

o

Apply a thresholding technique to binarize the results

Compute the area covered by the ltered patterns for each angle band Com

pare the results with the known areas of the directional patterns

Discuss the limitations of the methods and the artifacts in the results