Electrical and Computer Engineering

ECE Technical Reports

Purdue Libraries Year 

Automatic Gradient Threshold

Determination for Edge Detection Using
a Statistical Model A Description of the
Model and Comparison of Algorithms
Peter V. Henstock,∗ David M. Chelberg†

∗ Purdue University School of Electrical and Computer Engineering,

† Purdue University School of Electrical and Computer Engineering,
This paper is posted at Purdue e-Pubs.
http://docs.lib.purdue.edu/ecetr/95
AUTOMATIC GRADIENT THRESHOLD
DETERMINATION FOR EDGE
DETECTION USING A STATISTICAL
MODEL: A DESCRIPTION OF THE
MODEL AND COMPARISON OF
ALGORITHMS

TR-ECE 96-13
AUGUST 1996
 Automatic Gradient Threshold Determination
for
Edge Detection Using a Statistical Model
A Description of the Model and Comparison of Algorithms

Peter V. Henstock, David M. Chelberg

School of Electrical Engineering
Purdue University
West Lafayette, Indiana 47907-1285

This work was supported in part by NSF Grant No. IRI-90114L21.

 Contents

I Introduction 1

I1 Histogram Modeling 3
A Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
B Our Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

I11 Description of the Different Methods 6

A Global Descent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
B Two Step Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
B.l Alpha-Beta Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 8
B.2 Percentage Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
B.3 Threshold Determination . . . . . . . . . . . . . . . . . . . . . . . . . . 11

IV Experiments 12
A Data Preparation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
B Running the Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
C Testing Other Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
D Computed and Subjective Thresholds . . . . . . . . . . . . . . . . . . . . . . 15
E Robustness to Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

V Conclusion 23
 List of Tables

1 The mean and variance of the average distance (area) between the computed
and subjective thresholds, averaged over the 16 images. . . . . . . . . . . . . 18
2 A list of the algorithms presented in Table 1 . . . . . . . . . . . ,. . . . . . . 18
3 The mean and variance of the area between the subjective and computed
thresholds of the 5 different alpha-beta estimation methods. . . . . . . . . . 19
4 The mean and variance of the area between the subjective and computed
thresholds of the 4 different percentage estimation methods. . . . . . . . . . 19
5 The mean and variance of the area between the subjective and computed
thresholds of the 3 different threshold estimation methods. . . . . . . . . . . 20
 List of Figures

A bimodal histogram which can be easily thresholded. . . . . . . . . . . . . . . 5

A typical histogram for edge thresholding with the estimated densities of the
non-edge and edge pixels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
With some histograms, noise is required to spread the peak near zero gradient
units.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . , , . . . . . . 12
A battery of 16 images of varying scenes. Moving left to right starting at the
top, the images are: a house scene, an xray skull image, bethl, contact, crowd,
dilts, a synthetic dragon cartoon, flir images of a truck, man (gih), girl2, jo,
john, a satellite image of the earth, Madonna quantized to a few gray levels,
a text scan, and a still life (im16). . . . . . . . . . . . . . . . . . . . . . . . . 16
The thresholded images using the best algorithm. . . . . . . . . . . . . . . . 17
Images thresholded at 27, 39, 57, and 128 gradient units for subjective deci-
sions. The computed threshold is the second image. . . . . . . . . . . . . . . 19
The percent of non-edge pixels in each of the 16 images as determined by
the average of the 5 researchers (solid lines) and the average of the 5 best
algorithms (dotted lines). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
Synthetic images with added white Gaussian zero-mean noise of variances 1.0,
128.0, 1024.0, and 4096.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
Type I, Type 11, and overall misclassification errors for different noise levels
using the best algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
Histograms for the synthetic image under three different added noise variances. 26
 Abstract

In the enhancement/thresholding method of edge detection, the gradient values of pixels

exceeding a certain threshold are designated as edge pixels. However, selecting a threshold
has commonly been performed through ad hoc measures. This paper describes a method for
automatically selecting a threshold using a 5-parameter model. The model is based on the
weighted sum of two gamma density functions corresponding to edge and ncln-edge pixels. A
variety of statistical and fitting methods for finding the model parameters were evaluated by
comparing their computed thresholds to perceptual thresholds determined by subjects for
16 different images. The performance of the model was also analyzed under different noise
levels.

Index Terms: Automatic thresholding, edge detection, gradient, Sobel operator, statistical
classification.
 I Introduction
The goal of most object recognition systems is to identify the major features of an image
and match them against given models of known objects. Shape is commonly used for iden-
tification since it is one of the most distinguishing characteristics of an object. The shape
of an object is generally considered to be a high level property that can be computed from
a hierarchy of various lower order structures. At the lowest level, edgels or edge pixels rep-
resent the edges in the image that are grouped into line segments and then into various
curves, lines, rectangles, etc. Making accurate decisions about edgels at the lowest level will
facilitate more rapid and better performance at higher levels of the system. Consequently,
the edge detection process is of fundamental importance in object recognition. This paper
describes a model-based method for accurately and automaticall y determining a threshold
that separates edge pixels from non-edge pixels in intensity images.

Edges are defined as sharp intensity changes over a small area of an image. Since this general
definition does not lend itself to a specific mathematical formula, many edge detection algo-
rithms have been used. Edge detectors fall into three major classes: edge fitting operators,
zero crossings of second derivative operators, and enhancement /thresholdiing operators [I].
The first category attempts to orient a mathematical model to the edge using a best-fit ap-
proximation. The second method attempts to find the inflection points of the edge that are
the zero crossings of the second derivative operator, by definition. This paper focuses on the
third category of the enhancement/thresholding methods since they are the most common
in practice [2, 31.

Enhancement/thresholding methods apply differential or gradient-based operators to an im-

age. The operator is typically a 3 x 3 or 5 x 5 matrix which is convolved with. the image. The
Sobel or Prewitt methods use two operators which measure the local differences in contrast
in the horizontal and vertical directions and combine the results into gradient units. Areas
of high contrast have a high gradient value; areas of low contrast, a low value.

A threshold of the gradient values must be used to classify each pixel of an image as either
an edge or a non-edge pixel. Consequently, even in this class of operators, the presence of
an edge is still imprecise. Many differential operators have been studied i.n depth, yet the
determination of a threshold is still very difficult since it may depend upon. the application,
the source of the image, and the subjective perception of the viewer.

Our approach to identify the threshold dividing edge and non-edge pixels is based upon a
statistical model of the gradient values. By creating a histogram of the gradient values, we
can fit a statistical density function to both the edge pixel values and non-edge pixel values.
Once the densities have been found, the problem of finding the threshold is equivalent to a
statistical classification problem. This paper analyzes the effectiveness of several approaches
in finding the density functions and the corresponding threshold.

Studies in the evaluation of different edge detectors have contributed indirectly to work in
automatic threshold determination. In [I, 4, 51, several edge characteristics were identified
for comparing edge operators:

1. the number of correctly identified edge points

2. the number of incorrectly identified edge points (false positives)
3. the number of missed edge points (false negatives)
4. the number of multiple detections for a single edge, which is related to edge thickness
5. the distance between the identified edge and the true edge
6. edge continuity
After applying different operators to test images, the above characteristicis were quanti t a-
tively assessed and combined into a weighted error sum, allowing for comparison. In order
to compare the operators' performance, a threshold had to be selected. Venlratesh and Kitch
[:I.] used a least-failure-measure (lmf) threshold which minimized their error vector. While
qualitative measurement of an error sum has desirable characteristics, it is not clear that
the minimization of the error sum will, in fact, lead to an accurate determination of the
threshold. McLean and Jernigan [6] selected thresholds according to the "visual quality
of the resulting edge map, not according to any objective measure of the actual operator
performance". Most of the work in edge evaluation has focused on synthetic images with
known edges or compared the results of the operators against themselves, rather than against
an absolute standard, or a model-based criteria [5]. As a result, these methods cannot be
used to select an edge threshold, but serve to validate our use of both synthetic images and
subjective thresholds to evaluate the performance of the edge thresholding.

Previous work in gradient threshold selection includes the work of Haddon [7] and Zuniga
and Haralick [8]. After forming an estimate of the noise within the image, Haddon derived
a probability density of the edge strength from which he could select a threshold. With an
accurate measure of the noise, a global threshold could be selected to be independent of the
strength and number of edge pixels in the image. However, the approach assumes uniform
noise and only one image demonstrates its performance.

Zuniga and Haralkick [8] approached threshold selection as a Bayesian decision problem
involving two densities. Rather than fit densities to the edge and non-edge pixels, they
derived two conditional densities based on the gradient histogram by usin.g a facet model
in conjunction with a hypothesis test of the gradient values. For synthetic images, they
computed thresholds as a function of the noise that were better than subjective decisions,
although no real images were tested.

In Section 11, we survey automatic classification methods and describe our model in de-
tail. Section 111 introduces the different techniques of statistical estimation used to predict
the parameters of our model. Section IV describes and presents the results of our exper-
iments comparing the statistically calculated thresholds with the subjective perception of
edge thresholds, and evaluating the performance of the calculated t hresholdis under different
noise conditions.
I1 Histogram Modeling

A Background

The decision between edge and non-edge pixels based upon a gradient histogram is very
similar to the automatic segmentation/thresholding problem. Using grayscale or color as the
distinguishing criteria, pixels can be classified automatically. In determining the threshold
for different classes, these methods rely upon the modality, shape, or moments of the classes
in the histogram [9]. For the bilevel case in particular, these algorithms generally define the
location of the threshold to be at the minimum value of the valley between two peaks of the
histogram, which can be clearly seen in Figure 1. As Kapur, Sahoo, and Wong [lo] point
out, different algorithms try to improve the histograms by making the threshold more visible
through statistical or other enhancement techniques.

Due to the nature of the gradient thresholding problem, most of these bilevel classification
techniques cannot be used. Histograms of the gradient image normally do not have two
visible peaks and a valley as shown in Figure 2. This means that using the modality, shape,
or moments to classify the peaks will be futile since the two peaks and the valley separating
them cannot be identified.

Other classifications such as the p-tile [9] assume that the percentage of pixels in one class
is known a priori. The percentage of edge pixels is usually estimated to be 10 to 20% of the
overall pixels of an image [3, 11, 121. However since these percentages depend strongly upon
the images and the signal to noise ratio (SNR) present, such techniques will only provide a
rough estimate of the threshold and not an accurate measure.

B Our Model

Our approach employs a statistical classification based upon a 5-parameter model. As men-
tioned before, the gradient value reflects changes in grayscale over a local region of an image.
This change is generally the result of both edges (contrast changes) and noise. We assume
that we can statistically model the data by probability density functions. Given a 512 x 512
image with over 250,000 points, we have a large sample of data points which can be used
to find density functions. We model edge pixels and non-edge pixels by using two density
functions which sum to the original histogram.

The density functions which we use to model the edge and non-edge pixels are gamma density
functions. We assume that the edge and non-edge contrasts, as defined by the horizontal
 Gradient Squared Histogram
0.012

0
Gradient Squared Units

Figure I: A bimodal histogram which can be easily thresholdled.

Histogram of the Gradient Squared

0.014

Gradient Squared

Figure 2: A typical histogram for edge thresholding with the estimated densities of the
non-edge and edge pixels.
density functions. The second step, "percentage estimation", computes the ratio to edge and
non-edge pixels in the image. This section describes these estimation techniques in detail.

A Global Descent

For the global descent method, the goal is to find the combination of five parameters that
best fits the gradient histogram data. We are searching in a 5 dimensional space for a set
of parameters that most closely approximates the histogram. We use a ]?owell algorithm
[13], a multidimensional minimization based upon successive line minimizaltions that do not
require a gradient calculation. It utilizes conjugate directions in forming line minimizations
which converge quadratically to the minimum value of the function.

In order to use this search strategy, we had to define a "best fit" of the model to the data.
We investigated three distance measures for the entire experiment, each of ,which provides a
measure of how well a model with given parameters fits the actual data points. The distance
measures identify the difference in quantized curves x and y.

d l : the absolute distance error given by C i( yi - xi I

d2: the squared error given by Ci(yi - 2;)'

d3: the area between the curves using a trapezoidal approximation

While the measures are strongly correlated, they represent different accuracies and compu-
tation speeds, and vary with the amount of noise present. For the global descent method,
only the d l measure was used since the method is already quite computationally expensive.

B Two Step Estimation

Breaking the estimation process into two steps, the alpha-beta estimation anld the percentage
estimation are computed alternately until a convergence is reached. This: strategy is the
same as the EM (expectation-maximization) algorithm formulation [14, 151 and significantly
reduces the search space. In alpha-beta estimation, the two parameters For each density
function are estimated (four parameters total), given p;

a,k,p,k,a::,plkI P;
where pk is the relative probability of density0 verses densityl at the kth iteration. With
these estimates, a new pt parameter is formed in the percentage estimatio:n step
PO
k+l 1 p k ,k
07 0 7 19 1 '
pk

These two estimation steps are repeated until the estimated values converge or until no
progress is being made.

Since it is difficult to accurately estimate the parameters for two densities simultaneously
from their sum as in the first step, we divided all the points in the histogram into two non-
overlapping groups using the prior estimates for the density parameters and the po value.
Again using the EM algorithm formulation, we compute a ratio of the densities weighted by
po at each point of the histogram, defined as

where density0 and densityl are the reconstructions of the gamma densities given the esti-
mates of the a and p parameters from the most recent estimation. It should be noted that
this ratio function is actually not a density function, but only a ratio of densities with values
ranging from 0 to 1. From this ratio function we form two separate density estimates

for which new a and parameters can be more easily estimated. By truncating the histogram
where possible, as described later, the computation of this step can be reduced since it is
directly proportional to the length of the histogram.

B.l Alpha-Beta Estimation

We compared five different methods for estimating the a and ,l3 parameters of the two
densities. Two of these employed an iterative estimation and the remainder employed a
two-parameter Powell descent algorithm using different distance measurements.

1) Method of Moments This method generates point estimators by equating the first
two moments of the data with the first two moments of the parameter estimates given by
the following equations:
1 "
Datamomentl = - zi C
;=I

Estimatemomentl = a4

Solving these for a and ,B, we get n

2) Maximum Likelihood Estimation As the name implies, this approach is the formal
ML estimation using the derivative of the gamma density function which finlds the maximum
of the product of densities using the derivative's zero crossings:
n n G l x j(-1) e - X 0x i
I-J f (xi I a d ) =
i=l panr(a)n
Since the log is a monotonic function, we can take the log of both sides

l o g ( n f (xi / 8 ) ) = ( a - 1) z C xi
log xi - -- a n log P - n log I'(a).
P
Taking the derivative of this with respect to cr and ,b' and setting the resulting equations to
0, we obtain the MLE of a and P, respectively:
d d
-+ log(xi) - n log 4 - n- log[r(a)] = 0
da da
where

and

It should be noted that the d / d a of the log of gamma function requires significant computa-
tion with the integral. While it may be possible to perform a table look-up of precomputed
values for similar images, this is impractical over a wide set of images due to the high
precision required and large range of possible alpha values.

3-5) The Powell Descent Method This method is the same overall algorithm as the
global descent method described earlier in this section. The difference is that the Powell
method for alpha-beta estimation is searching for two parameters for each density, and each
density is handled separately. The search space is significantly reduced frorn the five dimen-
sions of the global descent method to only two dimensions. The Powell method minimizes
the distance between the model of the density and the density computed from the histogram
ratio function. That is,

We used three different Powell methods, for alpha-beta algorithms 3-5, with respective dis-
tance functions d l , d2, and d3 as defined in Subsection A.

B.2 Percentage Estimation

Four different algorithms were compared for the percentage estimation of non-edges. The
percentage of non-edge pixels or po is a scaling ratio parameter for the two densities. It
determines the relative weights of the two densities which sum to the gradient histogram.
Three of the algorithms use the Golden Section descent method with the distance measures
d l through d3. The fourth employs the EM algorithm [14, 151 on a Bernoulli estimate of
the ratio.

The Golden Section [13] descent algorithm is a one-dimensional minimization method based
upon bracketing the minimum value, that is quite efficient in searching for a minimum value.
Given the two densities described by the a and P parameters, po was defined as the value
which minimized the distance between the histogram and the sum of the two ;gammadensities
weighted by po. The equation is:

d i s t a n c e { h i s t o g r a m ( i ) , [densityO(i)x po + densityl(i) x ( 1 - po)])

i
where the distance is one of the aforementioned distance functions d l through d3.
The fourth estimation algorithm is based on a Bernoulli estimate of the ratio po and uses
the EM algorithm. The equation is:
ksum
Po = -
n
where

ksum = histo(i)
po xi densityO(j)
i +
po C j ( d e n s i t y O ( j ) ) ( 1 - po) Ci d e n s i t y l ( j )
and n is the total number of points in the histogram.

B.3 Threshold Determination

Given two overlapping densities, there are several methods of determining a threshold de-
scribed in pattern recognition and signal processing texts. Most of these methods are de-
signed for normal densities, rather than the gamma densities of our approach. Although
some researchers such as Canny prefer dual thresholds for edges, this paper will focus on sin-
gle threshold estimation. The results can be generalized for dual thresholds with our model
of the edges with an appropriate threshold criteria given the edge model. U'e considered the
MAP and ML decision thresholds, which are the most common in practice, and also tested
a direct po threshold. These three methods of computing the threshold were compared in

the experiments for their accuracy using our model.

The MAP threshold decides that an edge is present when

p ( z I e d g e ) x p(edge) > p(a: 1 non-edge) x ~(non-edge)

and a non-edge otherwise. In the notation of this section, we can rewrite this as
d e n s i t y l x ( 1 - po) > d e n s i t y 0 x po.
The maximum likelihood threshold decides that an edge occurs where
p(a: I e d g e ) > p(a: ( non- edge) or densityl > densityo.
This differs from the MAP decision in that it ignores the prior probabilities of edge and
non-edges.

The p-tile method of selecting edges forms a threshold at a given percentage, ignoring all
other information. Since po is a parameter of the model, we tested this approach, even
though this threshold is based directly on one parameter rather than the 3 parameters of
the MAP threshold and the 2 parameters of the maximum likelihood threshold.

IV Experiments

Two different techniques were used to analyze the effectiveness of our model. In the first
experiment, 16 diverse images were processed with the different algorithms and the computed
thresholds were compared with the subjective edge threshold decisions made by 5 researchers.
The second experimented evaluates the robustness of the model to different noise levels
applied to a synthetic image. We will discuss the image preparation and the procedures in
computing the thresholds, followed by the two evaluation methods.

A Data Preparation
All the methods described in Section I11 operate on the gradient histogram prepared from
the images. To standardize the testing, all the images were prepared using the following
steps: adding noise, smoothing, applying the Sobel operator, and truncating the histogram
which will be explained in more detail.

Histwram of Gradient Squared

Figure 3: With some histograms, noise is required to spread the peak near zero gradient
units.

Additive white Gaussian noise N(O, 1) (zero mean and unit variance) was added to each
image and rounded to the 255 grayscale values of the original images to spread the sharp
spikes which occurred in some of the histograms. For images with very low noise and little
texture which often occur in synthetic images, the histogram takes the form of a delta
function at gradient value 0 (no edge), and a small peak at a gradient valu'e proportional to
the contrast of the edge shown in Figure 3. Recalling that the gamma density at x = 0 is
zero, it is clear that the pixels of the non-edge density will not accurately fit the data. By
adding Gaussian noise N(0,l) to the image, the gradient values are shifted away from the
zero value, spreading the delta function. Since the resultant thresholds are robust to the
addition of low variance noise, this noise was added to all the images.

The image was then smoothed using a 3 x 3 Gaussian filter with unit varia.nce. A Gaussian
filter is commonly applied before edge detection to reduce the effects of noise by smoothing,
and as a form of regularization as described by Poggio [16]. The smoothing had virtually no
effect on the overall shape of the histogram.

Although our algorithm can be used by any differential or gradient based operator obtaining
similar results, only the 3 x 3 Sobel operator was used in this experiment. This operator was
chosen since it is effective in the presence of noise and is widely used [3, 21. After convolving
the Sobel operator with the image, the output gradient values were placed in a histogram
with a bin size of 1 gradient unit.

The largest 0.01% of the gradient values of the histogram were truncated which served several
purposes. First, by reducing the length of the histogram, computation wa;s saved since the
I

statistics describing the data fit need only be computed for pertinent values. Second, by
eliminating these outlying values, the accuracy was improved since such values tended to
bias both the statistical-based met hods and the fitting methods. Using different percentages
for truncation did not change the resulting thresholds as long as the percentages were much
less than the percentage of edge pixels in the image. The histogram was then normalized to
unit area for purposes of comparing data sets.
B Running the Algorithms
The goal of this portion of the experiment was to find and compare thresholds from the differ-
ent algorithms under the constraints of our model. We applied each method to the prepared
histograms described earlier in this section. An initial seed of 88% non-edge pixels was used
for all the images since they commonly compromise between 85% and 90% of image pixels,
though this is strongly dependent upon the image and the definition of edges. For all cases
except the global Powell descent involving simultaneous estimation of all the parameters, the
first iteration of the parameter estimation step utilized the maximum likelihood estimation
(MLE) method. This step provided a standard, quick, and reasonable initialization using
the updated percentage information for all the algorithms and did not appear to bias any of
the algorithms.

Following the initial estimation step, subsequent parameter estimation followed by the per-
centage estimation was performed and the process repeated as shown:

1. initial p: = 88%
2. initial alpha-beta estimation (p:) returns a,",p,", a!, P:
3. overall estimation loop (
percentage estimation (a!, pi, a!, Pf")updates JI;+'
alpha-beta estimation (a:, ph, a:, P:, p;") returns a!", p;+', a!+', P:+'
k=k+l)
where the k denotes the iteration number.

The complete iterative process of parameter estimation and percentage e~~timation,

shown
in the overall estimation loop was repeated a maximum of 12 times or until the threshold
converged within 6 gradient units of the previous iteration for all the algorikhms. The latter
condition required the minimum precision of the threshold to be approximately 0.5% (6 out
of a maximum 1200 gradient units) between on successive iterations.

The parameter estimation algorithms iterated for a maximum of 50 iterations or terminated

when the both the and PI parameters converged to within .3% of the pi:evious iteration.
The a parameters were not used as a termination condition since they were less sensitive
than the p parameters and also tended to vary proportionally to the ,B parameters. The
percentage estimation algorithms terminated when consecutive estimates (converged within
0.3%.

For all the methods, the range of valid percentage of non-edge pixels was restricted to
[1%,99%]of the truncated histogram. According to the 5-parameter model, the histogram
consisted of the sum of two non-zero gamma densities representing the ed.ge and non-edge
pixels. Bounding the size of the non-edge histogram away from 0% and 100% barred the
elimination of one of the edge densities.

C Testing O t h e r M e t h o d s
Before analyzing the results of our 5-parameter model, we tested several other techniques
described by Sahoo, Soltani, and Wong in their overview of thresholding methods, to verify
that these methods would not, in general, be effective for different gradient images. We
prepared the data using the same process of adding noise, filtering, and applying the Sobel
operator as before. As expected, the p-tile, node, and concavity methods vvere not effective
due to the varying shapes and edge to non-edge ratios of the different image histograms.
The p-tile method always selected the same percentage of pixels as the threshold, but the
percentage of edges in the battery of 16 images varied widely. The centroid and other
methods including the Otsu method could not find a threshold with the absence of two
peaks in most of the gradient histograms.

D C o m p u t e d a n d Subjective Thresholds
The purpose of the subjective analysis was to find how well the mathennatical threshold
compared with the subjective perceptual determination of edges. We used tihe human visual
system as the standard of edge detection since humans are able to rapidly and accurately
perceive global structures in images such as edges through a complex process, often referred
to as perceptual organization [17, 18, 19, 201.

For the experiment, the 16 images shown in Figure 4 were prepared as dlescribed in Sub-
Figure 4: A battery of 16 ima~gesof varying scEenes. Moving left to right starting at the top,
the images are: a house scene,, an xray skull i mage, bethl, contact, crowd, dilts, a synthetic
dragon cartoon, Air images a~f a truck, man (gih), girl2, jo, john, a satellite image of the
earth, Madonna quantized to a few gray levelIs, a text scan, and a still life (im16).
section A to form a histogram. Since the full range of thresholds is not necessary, we first
restricted the valid range of thresholds and then divided it into 256 equa1l:y spaced possible
thresholds, for greater accuracy.

Five subjects with experience in edge detection and computer vision were shown two black
and white 5.4" x 5.4" images on the computer screen-the original image and the thresholded
image. For each of the 16 images, the subjects were asked to raise or lower the threshold
using the mouse buttons so that the thresholded image matched their idea of a "best" edge
for the original image.

Figure 5: The thresholded images using the best algorithm.

Each of the algorithms of Section I11 were run on the histograms, giving a set of computed
thresholds. We used the area under the normalized histogram between each computed thresh-
old and the average subjected threshold as the metric for comparing the performance of the
different algorithms. Table 1 shows the mean and variance of the area differences for the
 Mean Var. a-,D % Est. Thr. Mean Var. a - ,8 % Est. Thr.
0.0611 0.0020 4 1 2 0.1381 0.0064 1 3 1
0.0618 0.0019 2 4 2 0.1387 0.0064 1 1 1
0.0633 0.0010 4 4 2 0.1387 0.0068 2 1 1
0.0741 , 0.0037 4 3 2 0.1399 0.0067 2 3 1
0.0836 0.0020 1 4 2 0.1520 0.0067 1 3 2
0.0963 0.0036 2 1 2 0.1522 0.0067 1 1 2
0.0967 0.0038 2 4 3 0.1522 0.0076 2 2 2
0.1016 0.0051 2 3 2 0.1562 0.0167 4 3 3
0.1084 0.0044 4 4 3 0.1577 0.0131 3 1 1
0.1166 0.0061 4 1 1 0.1577 0.0131 3 3 1
0.1219 0.0024 1 4 1 0.1606 0.0081 5 2 1
0.1229 0.0045 5 2 2 0.16180.0368 2 ;3 3
0.1234 0.0161 4 1 3 0.1620 0.0370 2 1 3
0.1324 0.0063
0.1362 0.0066
0.1365 0.0074
4
2
4
4
4
3
1 :I 1
1 1 1
0.1621 0.0042
0.1637 0.0124
1
3
:2
2 1
2
1
0.1651 0.0073 global glolbal global
1
Table 1: The mean and variance of the average distance (area) between the computed and
subjective thresholds, averaged over the 16 images.

I
Al~ha-BetaEstimation
I
11
I1
I Percentage
1 " Estimation
I
a-O I % Est.
1 Method of Moments 1 Golden Section (dl)
2 Max Likelihood Est. 2 Golden Section (d2)
3 Powell Descent (dl) 3 Golden Section (d3) 3 p-tile
4 Powell Descent (d2) 4 Bernoulli EM
5 Powell Descent (d3)

Table 2: A list of the algorithms presented in Table I

Figure 6: Images thresholded at 27, 39, 57, and 128 gradient units for subjective decisions.
The computed threshold is the second image.

best algorithms. Converting the thresholds into percentages of non-edge pixels, the average
of the 5 best algorithms and average subjective edge decisions are graphed in Figure 7. As
an example of the edges, three of the researchers chose the thresholds of 27, 57, and 128
for the im16 image shown in Figure 6. The best algorithm selected the second image of the
sequence with a threshold of 39. The thresholds for all the images with the best algorithm
is shown in Figure 5.
Alpha-Beta Estimation Method Mean Variance
1 Method of Moments 0.210 0.025
2 Maximum Likelihood Est. 0.142 0.013
3 Powell Descent with d l 0.257 0.026
4 Powell Descent with d2 0.262 0.084
5 Powell Descent with d3 0.313 0.042
Table 3: The mean and variance of the area between the subjective and computed thresholds
of the 5 different alpha-beta estimation methods.
1 Percent Estimation Method I Mean I Variance I

Table 4: The mean and variance of the area between the subjective and computed thresholds
of the 4 different percent age estimation methods.

For comparing the algorithms, Table 1 demonstrates that many of the a1go:rithms have very
similar performance over the different images. To better view the data, we have grouped
the data by averaging the values for every alpha-beta estimation, percent estimation, and
threshold technique in Tables 3, 4, and 5, respectively.
 "0 2 4 6 8 10 12 14 16
Image Number

Figure 7: The percent of non-edge pixels in each of the 16 images as deltermined by the
average of the 5 researchers (solid lines) and the average of the 5 best algorithms (dotted
lines).

I Threshold I Mean I Variance I

1 1 1 MLE
I I

1 0.221 (
I

0.037
1 2 1 MAP ( 0.187 1 0.029 1
3 1 p-tile 1 0.302 1 0.051 1
Table 5: The mean and variance of the area between the subjective and computed thresholds
of the 3 different threshold estimation methods.
The first two alpha-beta algorithms of Table 3 employ statistical methods to determine the
parameters, but the last three use the Powell Descent fitting technique. Although the results
in terms of accuracy and robustness are fairly similar, the statistical methodls are many times
faster than their counterparts. Consequently, the maximum likelihood estiination technique
is the preferred method in terms of both speed and accuracy.

In Table 4, the Bernoulli EM algorithm was the most accurate method for percent estimation
and was also much faster than the Golden Section counterparts. The reasoning is similar to
the statistical versus the fitting methods of the alpha-beta algorithms.

For Table 5, the percentage estimation step for the thresholding algorithms was not very
accurate as expected, since it does not use all the model information to fi.nd an optimized
threshold, but rather only one parameter of the model. The MAP thresholtl proved the best
method, but the MLE was also effective.

E Robustness to Noise

The previous section demonstrated that the 5-parameter model and corresponding alge
rithms were robust over a wide range of images. This experiment tests the robustness of the
algorithms over different noise levels applied to a synthetic image. Using the top algorithms
from Table 1, we added N(0, a ) noise to ten copies of a synthetic image with a 2 ranging from
1.0 to 4096.0, in multiples of 2.0. Four examples of the noisy images are shown in Figure 8.
We then calculated the threshold for each image at each noise level.

To analyze the data, we wanted to calculate the misclassification errors associated with the
average threshold choices at each noise level. Since the gradient values increiase with the noise
level (Figure 10)) the thresholds must also increase accordingly for the model to work. For
this evaluation, we assumed the a 2 = 1.0 case corresponded to the perfect c;lassification. We
then calculated the Type I, Type 11, and total errors as a function of the noise variance which
are plotted in Figures 9, respectively. As can be seen from the graph, the total classification
error is low (< .l) until u2 = 512 which corresponds to a SNR as low as 3.0. The reason
Figure 8: Synthetic images with added white Gaussian zero-mean noise of variances 1.0,
128.0, 1024.0, and 4096.0.

Noise Variance

Figure 9: Type I, Type 11, and overall misclassification errors for different noise levels using
the best algorithm.
that the total error is less than the Type I1 error is that the total error represents a weighted
sum of the errors. The relatively small error over the low SNR indicates that the algorithm
is effective at identifying a threshold. At high noise levels, the 3 x 3 Sobel is much less
effective as an edge detection operator which degrades the performance of any thresholding
technique. As an interesting note, the Type I1 error was significantly higher than the Type
I error. This implies that the researchers preferred more false positives tha~nfalse negatives
in estimating their thresholds.

V Conclusion

In object recognition systems, choosing a threshold to accurately find edge pixels in an image
at the lowest level can lead to significant computational savings at higher levels. Since most
automatic thresholding techniques do not apply to the specific problem of edge detection,
heuristic approaches are commonly used in research. We have developed ;a model of edges
that strongly agrees with the subjective perception of edges consisting of the weighted sum
of two gamma densities to represent edge and non-edge pixels. The modell proved effective
over a wide range of images and performed well in the presence of noise with a SNR as low
as 3.0.

After testing a number of the algorithms to calculate the thresholds based on the model,
we recommend both the Powell descent with d2 distance measure and the maximum likeli-
hood estimate in conjunction with the Bernoulli EM and MAP threshold. 'These algorithms
performed the best in terms of accuracy and computational speed for our :set of images.

This paper has discussed the use of automatic threshold determination over the entire image.
Many applications however require local rather than global thresholds. Our method can be
applied directly to such applications without modifying the algorithm by slimply taking the
gradient histogram data from selected regions of the image instead of the whole image. As
mentioned earlier, the model is also generalizable to operators other than ithe Sobel.

References

[l]S. Venkatesh and L. J. Kitchen, "Edge evaluation using necessary components," CVGIP:
 Graphical Models and Image Processing, vol. 54, pp. 23-30, 1992.
[2] B. K. P. Horn, Robot Vision. New York: McGraw-Hill Book Company, Inc., 1985.
[3] D. H. Ballard and C. M. Brown, Computer Vision. New Jersey: Prentice-Hall, Inc.,
1982.
[4] I. E. Abdou and W. K. Pratt, "Quantitative design and e ~ a l u a t ~ i oof n enhance-
ment/thresholding edge detectors," Proceedings of the IEEE, vol. 67, no. 5, pp. 753-763,
May 1979.
[5] L. Kitchen and A. Rosenfeld, "Edge evaluation using local edge coherence," IEEE Trans-
actions on Systems, Man, and Cybernetics, vol. SMC-11, no. 9, pp. 597-605, September
1981.
[6] G. F. McLean and M. E. Jernigan, "Hierarchial edge detection," C!omputer Vision,
Graphics, and Image Processing, vol. 44, pp. 350-366, 1988.
[7] J . F. Haddon, "Generalised threshold selection for edge detection," Putt ern Recognit ion,
vol. 21, no. 3, pp. 195-203, 1988.
[8] 0. A. Zuniga and R. M. Haralick, LLGradientthreshold selection using the facet model,"
Pattern Recognition, vol. 21, no. 5, pp. 493-503, 1988.
[9] P. K. Sahoo, S. Soltani, and A. K. C. Wong, LLAsurvey of thresholdling techniques,"
Computer Vision, Graphics, and Image Processing, vol. 41, pp. 233-260, 1988.
[lo] J . N. Kapur, P. K. Sahoo, and A. K. C. Wong, "A new method for gray-level picture
thresholding using the entropy of the histogram," Computer Vision, Graphics, and
Image Processing, vol. 29, pp. 273-285, 1985.
[ll] L, S. Davis, "A survey of edge detection techniques," Computer Gralphics and Image
Processing, vol. 4, pp. 248-270, 1975.
[12] M. Tanaka and T . Katayama, "Edge detection and restoration of noiisy images by the
expectation-maximization algorithm," Signal Processing, vol. 67, no. 5, pp. 753-763,
May 1989.
[13] W. H. Press, S. A. Teukolsky, W. T . Vetterling, and B. P. Flannery, Arumerical Recipes
in C. Cambridge: Cambridge University Press, 1992.
[14] M. Tanaka and T . Katayama, LLEdgedetection and restoration of noisy images by t h e
expect atons-maximization algorithm," Signal Processing, vol. 17, pp. 213-226, 1989.
[15] A. P. Dempster, N. M. Laird, and D. B. Rubin, "Maximum likelihood for incomplete
d a t a via the EM al g orithm," J. Royal Stat. Society, vol. B-39, pp. 1-3'8, 1977.
[16] V. Torre and T . A. Poggio, "On edge detection," IEEE Transactions on Pattern Analysis
and Machine Intelligence, vol. PAMI-8, no. 2, pp. 147-163, March 1986.
[17] A. Sha'ashua and S. Ullman, "Structural saliency: The detection of globally salient
structures using a locally connected network," Proceedings of the Second International
Conference on Computer Vision, 1988, pp. 321-327.

1181 A. Rosenfeld, "Pyramid algorithms for efficient vision," in Vision: Coding and Eficiency
(C. Blakemore, ed.), pp. 423-430, Cambridge: Cambridge University Press, 1990.

1191 J. R. Pomerantz and M. Kubovy, "Theoretical approaches to perceptual organization:

Simplicity and likelihood principles," in Handbook of Perception ana! Human Perfor-
mance (K. R. Boff, L. Kaufman, and J. P. Thomas, eds.), vol. 2, pp. 36.1-36.46, New
York: Wiley, 1986.

[20] A. Rosenfeld, "Theoretical techniques: Pyramid algorithms for perc:ept ual organiza-
tion," Behavior Research Methods, Instruments, and Computers, vol. 1:3,no. 6, pp. 595-
600, 1986.
 Gradient Value

Figure 10: Histograms for the synthetic image under three different added noise variances.