Multiple Resolution Segmentation of Textured Images.
A multiple resolution algorithm is presented for segmenting images into regions with differing statistical behavior. In addition, an algorithm is developed for determining the number of statistically distinct regions in an image and estimating the parameters of those regions. Both algorithms use a causal Gaussian autoregressive model to describe the mean, variance, and spatial correlation of the image textures. Together, the algorithms can be used to perform unsupervised texture segmentation. The multiple resolution segmentation algorithm first segments images at coarse resolution and then progresses to finer resolutions until individual pixels are classified. This method results in accurate segmentations and requires significantly less computation than some previously known methods. The field containing the classification of each pixel in the image is modeled as a Markov random field. Segmentation at each resolution is then performed by maximizing the a posteriori probability of this field subject to the resolution constraint. At each resolution, the a posteriori probability is maximized by a deterministic greedy algorithm which iteratively chooses the classification of individual pixels or pixel blocks. The unsupervised parameter estimation algorithm determines both the number of textures and their parameters by minimizing a global criterion based on the AIC information criterion. Clusters corresponding to the individual textures are formed by alternately estimating the cluster parameters and repartitioning the data into those clusters. Concurrently, the number of distinct textures is estimated by combining clusters until a minimum of the criterion is reached.
Introduction
The objective of texture segmentation is to separate an image into regions of distinct statistical
behavior. An implicit assumption in this process is that the statistics of each region are stationary
and that each region extends over a significant area. Therefore, it is reasonable to assume that
image pixels which are spatially close are likely to be of the same texture. In addition, any texture
segmentation algorithm which independently classifies each pixel in an image is likely to perform
poorly since locally there may not be sufficient information to make a good decision. For these rea-
sons, most segmentation algorithms either implicitly or explicitly impose some form of smoothness
in the resulting segmentation. For example, this may be done by dividing the image into arbitrary
blocks and classifying the texture of each block separately. However, if the block size chosen is
too small, discriminating among similar textures may be difficult. Alternatively, if the block size is
too large, regions of differing texture may be lost. In either case, the resulting boundaries will not
be accurate since there is no reason to believe that the actual texture boundaries occurred along
the block boundaries. Alternatively, a number of authors have proposed more natural methods of
imposing smoothness constraints on the segmentation of an image [6, 4, 3, 2, 1, 11]. These methods
use a random field with smooth spatial behavior to model the discrete valued field containing the
classification of each pixel in the image. Segmentation is then approached as a statistical estimation
problem. In this approach, a region in the image is classified differently from the surrounding
regions if there is sufficient statistical evidence to justify a distinct region regardless of size.
These methods of segmentation use Markov random fields (MRF) to model the discrete field
containing the individual pixel classifications. In an MRF, each point is statistically dependent only
on its neighbors so that the complexity of the model is restricted. The image is then segmented
by finding an approximate maximum a posteriori (MAP) estimate of the unknown field since
this is the most likely segmentation given the observed data. However, the resulting functions are
quite difficult to maximize globally. Stochastic relaxation algorithms such as simulated annealing
have been applied to these problems [6, 2]. It has been shown that, with the proper annealing
schedule, SA will converge to the global optimum [6], but such a schedule converges much too slowly
to be practically useful. In practice, an annealing schedule is chosen to yield a reasonable trade
off between performance and computation, but SA remains a computationally intensive method of
minimization.
Alternatively, a greedy minimization algorithm known as estimation by iterated conditional
modes (ICM) [4] is computationally efficient. This method maximizes the probability of the segmentation
field by deterministically and iteratively changing pixel classifications. When textures
can be discriminated over small regions containing few pixels, ICM yields good results. However,
our work indicates that when larger numbers of pixels are necessary to discriminate among different
textures, as is likely in higher resolution images, ICM is prone to being trapped in local minimum
of the cost function.
The algorithm developed in this paper first segments the image at coarse resolution and proceeds
to progressively finer resolutions until individual pixels are classified. At each resolution a greedy
algorithm is used to perform the classification, and then the result is used as an initial condition at
the next finer resolution. It is found that this multiple resolution segmentation (MRS) algorithm is
less likely to be trapped in local minima than ICM since at each resolution the regions of appropriate
size are classified and used to guide finer resolutions. This is consistent with the observation
that the convergence of MRF based segmentation algorithms can be improved by varying the
local neighborhoods [3]. In addition, the MRS algorithm requires substantially less computation
than ICM since constraints can propagate rapidly across the image at coarse resolution. These
advantages can be important in high resolution images where individual pixels contain relatively
little information, and in the discrimination between textures where information may be spread
over large regions.
The approach developed in this paper is analogous to the multigrid methods used in computing
solutions to partial differential equations. Multigrid methods, which have recently been applied to
some computer vision problems [13], substantially reduce computation since constraints propagate
rapidly at coarse resolution. However, when solving linear partial differential equations, multigrid
methods do not improve the quality of the solution since local algorithms are guaranteed to converge
to the unique solution (if finite precision effects are ignored). In contrast, multiple resolution
segmentation techniques can both reduce computation and improve performance since local minima
in the cost function may be avoided by allowing coarse resolution solutions to guide finer solutions
[11, 12].
In order to use algorithm in an unsupervised fashion, a method is needed to estimate the
number of texture regions and their parameters. Our approach is to first develop a statistical model
of the observed image based on the texture model used in segmentation. This overall model is
parameterized by the number of textures, the parameters of the individual textures, and the rate at
which each texture occurs. It is then shown that, for a fixed number of textures, approximate MLE
parameter estimation may be performed by alternately reestimating the texture parameters and
repartioning the clusters associated with each texture. This operation is analogous to algorithms
such as the K-means algorithm [17] but maximizes the log likelihood instead of minimizing the
mean square deviation about cluster means.
The problem of estimating the number of textures is quite different from estimating their parameters
since the MLE of the number of textures may be arbitrarily large. It has been suggested
[15, 16] that a reasonable alternative to the MLE is to find the number and type of regions which
minimize an information criteria (AIC) proposed by Akaike [14]. The AIC contains an additional
term to account for the bias due to the number of parameters to be estimated. We adopt a similar
approach, but based on a modified criteria. We also show that this criteria may be minimized
through the selective agglomeration of texture clusters and the application of the parameter estimation
algorithm described above. In this scheme, the decision to agglomerate two clusters is based on
a cost function which equals the change in the criteria when nodes are merged. This cost function
is shown to be a function of the texture statistics associated with each of the two clusters. As a
final step, segmentation is performed at coarse resolution, and data occurring along boundaries are
removed so that superfluous textures formed at boundaries can be eliminated.
Both the segmentation and parameter estimation algorithms use a nonhomogeneous Gaussian
autoregressive (AR) model for the individual textures. The model allows the extraction of a texture
statistic for each pixel which can be averaged over arbitrary regions of pixels to yield an aggregate
statistic for the region, and is therefore well adapted to use with multiple resolution techniques.
The advantages of a predictive texture model over a model which uses only mean and variance is
shown to be particularly important in estimating the number of distinct textures when the textures
have high spatial correlation.
Statistical Model
Throughout this paper, we will assume that the observed image, Y , is a random field defined
on a rectangular grid, S, of N points, and the value of Y at a point s 2 S will be written as Y s .
When necessary the points in S will be explicitly written as integer pairs (i; j). X will denote the
field, also of N points, which contains the the classification of each pixel in Y . Points in X will
take values in the set f1; ::; Mg, where M is the number of textures or classifications. Further, the
conditional density function of Y given X, is assumed to exist and to be strictly positive and is
denoted by f(yjx), and the probability P
Using this framework, the image may be segmented by estimating the pixel classifications X
given the observed image Y . In particular, we will use the MAP estimate of X.
Since the set of possible segmentations is discrete, the maximum or minimum of a well defined
function on the set must exist. Hence, once the distributions for f(yjx) and p(x) are defined, the
problem of segmenting an image will reduce to that of minimizing a cost function.
2.1 Texture Model
Given a known segmentation, x, the image Y will be modeled as a causal nonhomogeneous
Gaussian AR random field [7]. This model has the advantage that the extraction of texture statistics
and the estimation of parameter values are direct and computationally efficient. The causal AR
model parameterizes the conditional distribution of a point in the field given past values of points
in a local neighborhood. Of course, the concept of a past is dependent on a specific ordering of the
image pixels (in this work raster ordering is used). A question which arises is whether the artificial
introduction of causality to an image will restrict the range of textures which may be modeled
accurately.
A partial answer to this question is found by considering other Gaussian random field models
such as simultaneous autoregressive (SAR) and Gauss Markov random field (GMRF) [8] which do
not require the introduction of a causal ordering. For Gaussian random fields, the relationship
among AR, SAR and GMRF models depends on the neighborhoods which are used. The AR model
is a special case of the SAR model obtained by constraining the prediction coefficients for future
points to be zero. Also, any SAR model may be equivalently described as a GMRF but with a
larger finite neighborhood. However, by choosing a sufficiently large neighborhood size, any of these
models may be used to approximate arbitrarily closely a Gaussian random field with well behaved
spectral density function [7]. Therefore, the major practical difference among these models is the
complexity required to achieve a sufficiently accurate model. Because SAR and GMRF models have
the disadvantage that maximum likelihood (ML) parameter estimation is difficult, we choose to use
the AR model. This difficulty in parameter estimation arises due to the complex dependence of the
distributions normalizing constant on the parameters [8] being estimated.
The image Y is assumed to consist of M textures each of which is described by a parameter
is the index of the particular texture. Each parameter vector, ' m , contains the
mean, -m , prediction variance, oe 2
m , and prediction coefficients, am , of that texture. The parameter
set used at the point s 2 S is then ' xs . We will assume raster ordering of the points in S, and for
to mean
denote element by element subtraction of the components of s and r.
A basic property of the AR model is that the prediction errors, formed by filtering with the
prediction coefficients, are independent. This property can be used to calculate the conditional
distribution of Y given x. The forward prediction error at s for a P th order AR model is given by
~
a xs (r) (Y
where there are P values of r for which a xs (r) may be nonzero. The placement of these P values
depends on the prediction window used. We will assume the use of a second quadrant predictor
window although other models could also be used. If we ignore the effects of the image boundary,
the components of the prediction error ~
Y s are independent with mean zero and variance
s
Given this information, the conditional distribution of ~
Y given x is just that of N independent
Gaussian random variables. However, the Jacobian of the transformation from Y to ~
Y is needed to
compute the conditional distribution of Y . This can be calculated by arranging the components of
Y and ~
Y into vectors with raster ordering. It is then possible to write (2) in the form
~
where A is a matrix containing the prediction coefficients and C is a constant vector. Also, A is
lower triangular, and its diagonal entries are all 1's. This is due to the fact that in raster ordering
all points used in prediction will precede the point being predicted. These facts imply that the
Jacobian of the transformation from Y to ~
Y is 1, and that the conditional distribution of Y is
Y
xs
exp
~
s
xs
This results in a negative conditional log likelihood of
~
s
xs
The log likelihood function (3) has the form of a sum of local functions of y and x. In general, the
methods used throughout this paper are applicable to any texture model which has a log likelihood
function of the form
l s (yjx s
where the functions l s depend on all of y but only x at the point s, and c(\Delta) is an arbitrary function
of y. Specifically, for this texture model and the local likelihood functions have the form
l s (yjx s
~
s
xs
2.2 Segmentation Field Model
We use an MRF to model X, due to its isotropic nature and its restriction to local interactions.
In order to define an MRF, the concept of neighborhoods must first be defined. The neighborhood
of a point s 2 S is a set of points @s ae S with the two properties that 8s; r 2 S s 62 @s and
@s. The set of ordered pairs f(s; @s)g s2S is then known as a neighborhood system.
We will use an eight-point neighborhood system in which the neighbors of interior points in S are
given by
The neighbors of boundary points will depend on whether a toroidal or free boundary is specified.
For an MRF, X, the conditional distribution of a point in the field given all other points is only
dependent on its neighbors; that is
clique, c, is a subset of points in S such that if s and r are two points contained in c then s and
r are neighbors. Notice that the set of all cliques is induced by the neighborhood system. For a
given neighborhood system, a Gibbs distribution is defined as any distribution, p(x), which can be
expressed in the form,
z
exp
all cliques c
where the functions V c are arbitrary functions of the values of x on the clique c, and z is a normalizing
constant. The constant T is physically analogous to temperature, and the argument of the
exponential
all cliques c
is physically analogous to energy.
The Hammersley-Clifford theorem [5] states that X is an MRF with strictly positive distribution
and a neighborhood system f(s; @s)g s2S if and only if the distribution of X can be written as a
Gibbs distribution with cliques induced by the neighborhood system f(s; @s)g s2S . Therefore, if p(x)
is formulated as a Gibbs distribution, we know that X will have the properties of an MRF. The only
cliques which will be used in the energy function for p(x) are pairs of horizontally, vertically and
diagonally adjacent points. The functions of these cliques will appear implicitly through the two
functions is the number of horizontally and vertically neighboring points
of different value in X, and t 2 (X) is the number of diagonally neighboring points of different value
in X. The density function for X is then assumed to be of the form
z
exp
are parameters of the model. Since both t 1 and t 2 may be written as the sum
of functions of cliques, we know that (6) is the distribution of an MRF with an eight-point neigh-
borhood. A more general model for MRF's with two point cliques results by allowing a different
potential for each distinct combination of unlike pixel pairs, and by adding a term for one pixel
cliques. It should be noted that the arguments used in this paper may be easily generalized for
this case. However, we will assume that prior knowledge is not available about the behavior of
individual class regions.
Substituting in p(x) from (6) and log f(yjx) from (4) into (1) results in the final form of the
minimization problem
where
l s (yjx s
3 Minimization Methods
In general, the problem of minimizing (7) is quite difficult since the cost function is not convex
and contains complex local minima. The problem may be formulated as a one dimensional dynamic
programming problem by associating a state with the value of x on a entire line, but the resulting
algorithm is computationally infeasible. Derin and Elliott have proposed finding an approximate
solution by solving the dynamic programming problem along thin strips [1]; however, the computation
required for such an algorithm increases dramatically with the number of textures. Also, while
the problem may be solved in one pass, this pass is quite complex and sequential in nature. The
approaches explored in this paper are iterative, but generally only require the application of simple
local operations, and they are well suited to highly parallel implementations.
3.1 SA and ICM Methods
SA is a general purpose algorithm for minimizing a cost function, C(x), as a function of the
state x. SA works by generating an irreducible Markov chain of these states. At each step of the
chain, the transition distribution is chosen so that the limiting distribution is a Gibbs distribution
with temperature T and energy function C(x). The minimum of C(x) is then found by allowing T
to drop slowly as the chain is generated. Since for sufficiently small T the lowest energy state of
the Gibbs distribution occurs with high probability, it is reasonable to expect that if T drops slowly
the samples will be concentrated about this minimum energy state. An advantage of SA is that it
only requires that the change in cost be computed between consecutive states in the Markov chain.
There are two distinct methods of generating samples from such a Gibbs distribution, and
consequently two distinct methods of performing SA. Both depend on the ergodic properties of
the Markov chains formed by iteratively replacing points in x. In the Metropolis algorithm [9] a
symmetric transition probability is defined and then transitions are either accepted or rejected with
a probability depending on the change in energy. The alternative method, which will be used in
this paper, is the Gibbs sampler [6] which replaces points in x with samples from the conditional
distribution of a point given the remainder of the field.
Notice that when C(x) is the energy function of an MRF then p T only depend on x @s .
For the segmentation problem formulated in the previous section the cost function is given by
U(xjy). Because the functions l s (y j\Delta) are only dependent on individual points in x, it is easily
verified that U(\Deltajy) is the energy function of a Gibbs distribution with an eight-point neighborhood.
Therefore, by applying the Hammersley-Clifford theorem, we know that the conditional distribution
of (8) will only be a function of the eight neighbors of x s . If T is large, transitions will occur almost
uniformly over the set of possible values for x s ; but if T is small, the value of x s corresponding to
the most likely conditional probability will be chosen with high probability.
It has been shown [6] that if T is varied according to a schedule, T (n), with the correct starting
value and rate of decrease, then as n !1 the distribution of X(n) will be uniform over the values
of x which minimize the cost function. In practice this annealing schedule is much too slow, so a
faster schedule is used which represents a trade-off between performance and computation. When a
faster schedule is used the algorithm is not guaranteed to escape from local minima. For example,
if a large region has been misclassified, the algorithm will probably never explore the possibility of
changing all of the classifications in that region. For this reason, the initial condition to the SA
algorithm is often of critical importance.
Besag has proposed the iterated conditional modes (ICM) algorithm [4] as an alternative to SA.
ICM may be regarded as SA with the extreme annealing schedule T replaces points
in x with a value which maximizes the conditional density function at each point. Since
arg min xs
it is clear that this replacement selects a value for x s which minimizes the cost under the constraint
of fixing the remaining values in x. Let v 1 be the number of horizontal and vertical neighbors
of x s which are not equal to k, and let v 2 be the number of diagonal neighbors of x s which
are not equal to k. The replacement can then be explicitly expressed as
arg min xs
The ICM algorithm only differs from SA with T when the argument which minimizes
is not unique. In this case, SA will replace x s with a value chosen with uniform probability
from the minimizing values. However, ICM selects the previous value of x s if it minimizes the
cost. Otherwise, ICM chooses pseudorandomly among the minimizing values. Consequently, each
time this replacement operation changes the value of a point in x, the cost, U(xjy), is assured to
decrease. Besag suggested these replacements be performed for a fixed number of iterations through
however, our approach is to continue replacing points in x until no additional changes occur.
Since U(\Deltajy) is strictly positive with a finite domain, the convergence of the algorithm in a finite
number of steps is assured. The minimizing state x which is obtained is a local minimum in the
sense that any other state which differs from x at only a single point must have energy greater than
or equal to U(xjy).
Because ICM can not perform back-tracking, the initial condition is critical. A reasonable choice
for an initial condition is the maximum likelihood estimate (MLE), or equivalently the estimate
obtained by dropping the cost terms due to the a priori probability. Because the functions l s (y j\Delta)
are only dependent on individual points in x, the MLE is easily written as
l s (yjk) (10)
The algorithm for ICM may now be stated explicitly as follows:
1. Calculate -
xMLE as an initial segmentation using (10).
2. Perform the local minimization defined by (9) at each pixel in a specified order.
3. If changes occur then repeat step 2, otherwise stop.
It should be noted that the specific order in which the replacements are performed is important.
If raster scan ordering is used, current pixel replacements are affected by the classification of the
adjacent pixels in the raster ordering. Since the prior distribution on X encourages pixels which
are adjacent to be classified similarly, this often has the effect of biasing the position of boundaries
in the direction of the raster ordering. It is not surprising that the solution obtained by ICM could
be affected by the update ordering since the minimum reached is only local and may depend on
details of the implementation. Another disadvantage of raster ordering is that it requires sequential
replacements, and therefore does not allow for a parallel implementation.
The solution to these problems is to break x up into groupings such that the members of each
grouping are conditionally independent of each other given the points in the remaining groups.
Because of this property, any updating done within a group does not depend on ordering, and all
points in a grouping may be updated in parallel. Such a division of x is known as a coding [5] and
depends on the neighborhood system used. An eight-point neighborhood results in four groupings.
So, it is possible to perform the updating in parallel using a maximum of N=4 processors.
An important aspect of a deterministic replacement algorithm such as ICM is that if the neighbors
of a point have not changed then it need not be updated. This scheme can be implemented
by keeping a flag at each point which indicates whether it is necessary to update the point. At
the start of the algorithm all flags are set. Then each time a pixel is updated its flag is cleared
and the neighboring flags are set if the pixel's value changed. Consequently, only the first iteration
requires the application of the replacement operation at each point, and at later iterations most
of the computation can be concentrated to regions where significant change is occurring. This can
substantially reduce computation.
3.2 Multiple Resolution Segmentation (MRS)
In this paper, we shall perform segmentation using a multiple resolution approach. The MRS
algorithm will begin segmentation at the coarsest resolution and use the result as an initial condition
for segmentation at the next finer level of resolution. Each resolution will correspond to a level in
a quad tree, so a lattice point at one resolution will correspond to four points at the next finer
resolution. This group of four points will be referred to as a block, and X (n) will denote the field
containing the classification of each of the blocks at resolution n. In general, quantities which vary
as a function of resolution will have superscripts in parentheses to denote the level, and a superscript
of 0 will denote the original quantity, so . The mapping from pixel coordinates at the n th
resolution level to points at the previous coarser (n th level can be formally written as
where b\Deltac denotes the greatest smaller integer function. Also, D k (s) will denote the composition of
the function D(\Delta) with itself k times. This structure is illustrated in Figure 1.
We may formulate local log likelihood functions for each resolution n by assuming that the fine
segmentation, x, is an interpolated version of x (n) . Specifically, if 8s 2 S, x
then applying
(4) results in
log f(yjx (n)
l (n)
s (yjx (n)
where
l (n)
l (n\Gamma1)
r
l (0)
and D \Gamma1 (r) is the set of the points at the next finer resolution which map to r. The following
section will discuss specific methods for computing these decimated likelihood functions.
At each resolution, our approach will be to solve an analogous segmentation problem to the one
formulated for the fixed resolution problem. The optimization criteria we will use for finding the
best segmentation at resolution n is
l (n)
where t (n)
count adjacent values of different type in x (n) , and fi (n)
1 and fi (n)
2 are the parameters
of the coarse resolution MRF. This criteria may be minimized by using the local replacement
operation
arg min
l (n)
The remaining question is how to choose these coarse resolution parameters. One possible
approach is to choose these parameters so that the cost functions at coarse and fine resolutions are
equal
[11]. In fact, this may be done by observing that if X (n\Gamma1)
This relation implies that the fine and course resolution cost functions will be equal if
resolutions n. By using these parameters, minimization of (12) corresponds to the minimization
of the original cost criteria (7) under the constraint that the solution be constant on the
appropriate sized blocks.
Coarse resolution minimization using the parameters given by (14) will effectively minimize
(7) if the correct segmentation is approximately block constant. However, this recursion for the
parameters has an undesirable property. It implies that the MRF models for coarser resolution
segmentations should have progressively higher spatial correlation, or alternatively, finer resolution
segmentations should have lower correlation. This, of course, runs counter to normal assumptions of
spatial coherence in images, and will tend to cause insufficient spatial correlation at finer resolutions
or excessive correlation at coarse resolutions.
A more reasonable approach is to assume that the spatial correlation is independent of the
resolution since it avoids the problem of excessive correlation at coarse resolutions. Also, this
assumption is appropriate when prior information is unavailable about the likely scale of regions in
the image. Therefore, in all experimentation, we will fix the parameters of the MRF as a function
of scale.
The L level MRS algorithm may be defined as follows:
1. Compute l (n)
(y j\Delta) for the levels using (11).
2. Compute the ML estimate of X (L\Gamma1) as an initial condition using (10). Set
3. Compute
x (n) by iteratively performing the local minimization of (13) until no more changes
occur.
4. If Otherwise, compute an initial value for - x (n\Gamma1) by replicating the values of -
x (n)
for each block in -
x (n\Gamma1) . Set return to step 3.
3.3 Decimation of Likelihood Functions
A question now arises of how the decimated likelihood functions, l (n)
s (y j\Delta), should be calculated.
The direct method is to explicitly calculate the functions for each texture, k, and sum those functions
over the proper blocks. However, if the number of textures is large, this may require considerable
computation and storage. Also, it presumes that the parameters of the textures are known. An
alternative method is to extract statistics for each block so that the functions l (n)
s (yjk) may be
computed both as a function of k and the model parameters ' k . The parameters at different
resolutions may then be calculated by decimation of those statistics from fine to coarse resolution.
Let A ae S be a region with NA points over which the likelihood function is to be summed. Then
the regions likelihood function, l A (yjk), may be expanded using (4) and (5).
~
s
~
The parameters of the k th texture model are '
We may write the log likelihood as
an explicit function of these parameters by defining three vectors. Let a k be a vector of prediction
coefficients in raster order. Let y s be the vector of corresponding points in y, and let 1 be a vector
of equal dimension with components of 1. For example, if with a quarter plane prediction
window the vectors have the following forms.
a
The autocorrelation, RA , and mean, mA , statistics for the region A can be defined using these
vectors.
The prediction error may now be written as a function of these statistics.
~
a t
a t
a t
By substituting this formula for the prediction error into the expression for l A (yjk), we may verify
that RA and mA are sufficient statistics.
l A
a t
If we define a statistic vector - may define a new function, L(\Delta; \Delta), such
that
Thus, if one can find a recursion for calculating the region statistics, - (n)
s , from fine to coarse
resolution, then the corresponding likelihood functions, l (n)
s (yjk), may be calculated from these
statistics. Let A and B be two disjoint regions, and let may be shown using
(15) that the statistics for C, - C , may be calculated from - A and - B . The commutative operation
of combining region statistics will be denoted by the symbol \Phi. Therefore
where the components of - C are
The statistics at coarse resolutions may now be calculated from those at fine resolutions using this
operation. If we use a notation for the repeated combining of region statistics, then an equation
similar to (11) may be defined.
This operation of combining region statistics will also be important
in the agglomerative clustering used in unsupervised parameter estimation.
Finally, we will require a method for calculating the MLE estimate of the parameter vector '
for a region with statistic -. In order to simplify these calculations we will approximate the MLE
estimate of - by the average of the components of m. Therefore, -
is a unit vector
with equal components. If the autocorrelation is decomposed into the following
then the MLE values of the components of ' are given by
Unsupervised Parameter Estimation
In order to estimate the number and type of textures without supervision, an algorithm is needed
which is independent of free parameters that may vary from image to image. The approach taken is
to base the estimate on a statistical model of how much variation is likely among different texture
samples. This results in an algorithm which is somewhat analogous to the well known K-means
algorithm [17], but incorporates an additional step for determining the number of distinct textures.
4.1 Clustering Model
The first step in developing this clustering method is to calculate the joint probability distributions
for the image and the pixel classifications given the texture parameters and the number of
textures. The Gibbs distribution used in the previous section to model the pixel classifications can
not be easily applied to the problem of unsupervised parameter estimation. The difficulty arises because
the unsupervised estimation problem requires that the distribution of X be parameterized by
the number of textures, M . However, the dependence on M is quite complex since the normalizing
constant of the Gibbs distribution is a complex function of M . In fact, normalizing constants for
Gibbs distributions have only been calculated in a few special cases (e.g. homogeneous Gaussian,
Ising). This intractability introduces a significant problem since the normalizing constant is critical
when comparing the likelihood of different clustering hypotheses. However, we would still like to
use our prior knowledge about the spatial correlation of pixel classes in determining the number
and types of textures present.
The approach we take is to first separate the image into blocks of equal size (normally at the
coarsest resolution used in the MRS algorithm). The prior knowledge about the spatial correlation of
pixel classes is then incorporated by assuming that each block contains only a single texture. Later,
we will remove the effect of blocks which fall on boundaries and may contain a mixture of textures.
In segmentation, we assumed that the classes of blocks had a Gibbs distribution in the form of (6)
with properly chosen fi parameters. The Gibbs distribution embodied the prior knowledge that
adjacent blocks are likely to have the same texture. In this section, we ignore information regarding
the spatial ordering of blocks. Conceptually, this ordering information may be ignored by assuming
that the blocks have been randomly reordered. This is not an approximation since we may choose
to ignore this information. However, our ability to distinguish similar textures will be diminished.
The advantage of this approach is that the weakest possible assumptions are made about the prior
distribution of the data, so that the clustering model may more accurately describe the observed
data.
When the order of image blocks is ignored, the only relevant (sufficient) statistics are the texture
statistics for each block, and the number of blocks of each texture. For fixed M , the Gibbs distribution
determines the distribution of the number of blocks of each texture. However, for a general
choice of M and fi, a simple form for this distribution is not known. Therefore, we approximate
this distribution using a multinomial distribution with unknown parameters. It should be emphasized
that this is not equivalent to assuming that the ordered block classifications are independent
random variables. In fact, block classifications may have high spatial correlation and still result in
a multinomial distribution for the number of blocks of each texture.
All clustering will be done using texture statistics extracted from image blocks normally at the
coarsest resolution level, L \Gamma 1. Let A k ae S (L\Gamma1) be the set of blocks at this level which have a
texture of type k for . Denote the number of elements in the set A k as W
and the number of elements in S (L\Gamma1) as j. Notice that A k will be random since it is a
function of the random field X (L\Gamma1) . W k is in turn a random variable since it is a function of A k . The
sets, fA k g M
partition the field X (L\Gamma1) into groups each of which has uniform texture. Therefore,
the random sets fA k g M
k=1 contain the same information as the random field X (L\Gamma1) , and the two
expressions can generally be interchanged in probability distributions. However, we have chosen to
use the partitioning sets since they are suggestive of the clustering techniques to be developed. The
additional constraints that
make the first M \Gamma 1 values of A k
and W k an unambiguous specification of the M th value. Denote as the vector of
texture parameters.
In appendix A, the joint distribution of Y and the unordered block classes represented by
k=1 is calculated. As described above, the the number of blocks of each texture, W k , is
modeled by a multinomial distribution with unknown parameters, ae k , and information about the
relative position of blocks in the lattice is ignored. The distribution is then given by
Y
where
the functions L(\Delta; \Delta) are defined in (17) and the abbreviated notation - s is used in place of - (L\Gamma1)
s to
denote the statistic for block s. It is also shown in appendix A that by summing over all possible
partitioning sets the distribution for y is given by
Y
which is recognizable as the distribution of W independent and identically distributed random
variables with a marginal distribution formed by a mixture of M distributions.
While we now have an accurate closed form distribution as a function of the parameters ' and
ae, the problem of estimating the number of parameters remains poorly defined. For example, if the
MLE texture parameters are estimated for the cases of both M and M then the fit of
the M+1 texture case will always be as good or better since it includes the case of M textures. For
this reason, ML estimation will generally choose a large or infinite number of textures to achieve the
best fit to the observed data. However, this solution is unacceptable since segmentation is intended
to partition the image into the minimum number of necessary distinct textures.
Zhang and Modestino [15, 16] have suggested that a reasonable alternative to the MLE is to
find the number and type of regions which minimize the AIC information criteria proposed by
Akaike[14]. The AIC contains an additional term to account for the bias due to the number of
parameters to be estimated.
The condition of identifiability, which necessitates that the Fisher information matrix is invertible
[18], is important in this application. Setting ae would seem to result in a distribution with one
less parameter. However, in this case the mean, variance and prediction coefficients corresponding
to texture M are no longer identifiable. Therefore, the number of identifiable parameters is reduced
by (where P is defined in section 2.1 as the number of AR prediction coefficients). Applying
this definition of AIC to (21) yields
This optimization criteria has a number of disadvantages. Calculation of (22) requires the use
of W statistics, one for each block. This contrasts with the calculation of (20) which only requires
the use of M statistics, one for each texture. Minimization of (22) for fixed M is difficult due to
the lack of known partitioning sets. Approximate minimization may be performed using the EM
algorithm [19, 20] or equivalently Baum's algorithm [21]. However, the minimization process is still
substantially complicated by the lack of class information, and does not result in an algorithm with
a simple agglomerative clustering interpretation.
For these reasons, we propose a modified criteria which estimates the missing class information
together with the number of textures and their parameters. This criteria results from replacing the
log likelihood of y given in (21) with the joint log likelihood of y and the partitioning sets given in
(20). This modified clustering criteria is
The objective of all cluster techniques developed in the following sections will be to minimize this
clustering criteria as a function of the parameters, f'; aeg, the partitioning sets, fA k
k=1 , and the
number of textures, M .
Finally, we make a number of observations regarding properties of the clustering criteria given
in (23). The manipulation of this criteria only requires the use of M statistics, and it will be seen
that the resulting minimization techniques are analogous to traditional agglomerative clustering
methods. By fixing the number of textures and the partitioning sets, minimization of this criteria
yields the ML estimates of the parameters. By fixing the parameters and the number of textures,
minimization yields the MAP estimate of the partitioning sets. The minimization of (23) is, of
course, only unique up to a permutation of the texture parameters and the corresponding partitioning
sets. We can make the solution unique by imposing an ordering on the parameter sets such
as ae k - ae k+1 . (If equality holds, then a finer ordering may be imposed using the ' k parameters.)
For a fixed M , the solution of this minimization problem may fall on a boundary of the parameter
space corresponding to ae This corresponds to the situation in which texture M never occurs
and no blocks are classified with that texture. In this case, the optimum number of textures will
be less than M since minimization for M 0 - M will also result in ae M
4.2 Estimation
The number of distinct textures and their parameters will be estimated using three procedures
each of which strictly reduces the clustering criteria of (23). Each of the procedures may be thought
of as minimizing with respect to one of M , ('; ae) or fA k g M
while fixing the remaining two quan-
tities. The first procedure reestimates the parameters ' and ae from the region statistics - A k . The
second procedure repartitions the set S by choosing the sets fA k g M
which minimize (23). These
two operations are analogous to the centroid estimation and data repartitioning in the K-means
algorithm. In this context, each of the sets A k contains the members of a cluster corresponding to
texture k. The last procedure is an agglomerative hierarchical clustering algorithm which reduces
M by combining pairs of sets A k when this minimizes (23).
Because the terms involving ' and ae can be separated in the log likelihood, the clustering criteria
may be minimized by each variable independently. Minimizing with respect to ' is equivalent to
calculating the ML estimate of ' k with respect to - A k for each texture k. The algorithm for doing
this is given in (19). The ML estimate of ' k will be denoted by -
'(- A k ). The ML estimate of ae may
be found by minimizing (23) with respect to faeg M
k=1 subject to the constraint that
The result is
Choosing partitioning sets A k is equivalent to choosing the field x (L\Gamma1) which minimizes (23) for
a fixed set of parameters. This may be done by choosing the class of each block to be the most
likely texture given the observed data.
Therefore, the best partition of the blocks is the partition which assigns each block to its most likely
class.
A third mechanism is required to minimize (23) with respect to the number of clusters, M . The
approach we take is to use agglomerative hierarchical clustering. Let c(k; l) be the change in (23)
caused by the merging of two clusters, A k and A l , into a single cluster, l . We again
assume that the remaining parameters, ('; ae) and fA are held constant. However, in
order to define c(k; l) the parameters of the new cluster, C, must also be specified. We choose the
ML estimate corresponding to C since this is the parameter vector which yields the minimum value
of (23) after merging clusters. By using this definition of c(k; l), we obtain the following expression
'(- C ) is the ML estimate of ' given the statistics - l for the new cluster C. By
substituting the MLE into (16), it is easily shown that for any statistic, - A ,
A is the ML estimate of the prediction error given the statistic - A , and NA is the number of
pixels in A. In the complete clustering algorithm described at the end of this section, the functions
c(k; l) are only computed when the parameter values for clusters A k and A l take on their MLE
values. Therefore, we may apply (25) to yield the following simplified expression for (24).
+NA l log
A l
From this formula, we can see that the only information required to calculate c(k; l) is the statistic
vector for the regions, - A k and - A l .
Since the global objective is to minimize (23), we are justified in merging any two clusters for
which c(k; l) is negative. This merging operation can be repeated until c(k; l) is greater then or
equal to zero for all clusters, A k and A l . In practice, it is often the case that many combinations
of clusters have negative c(k; l). Since the minimum found by this agglomeration procedure is not
global, the order in which clusters are merged can affect the final result.
A steepest descent strategy first merges the two clusters with the minimum value of c(k; l).
A brute force approach to steepest descent requires M 2 calculations of c(k; l) to find the two
nearest clusters. However, the computational complexity of steepest descent may be reduced by
using balanced trees [22]. For each cluster, A k , a balanced tree is used to maintain a list of all
other clusters in order of increasing c(k; l). For each tree, insertions, deletions and searches may
be performed in O(log(M \Gamma 1)) time (since there are nodes in each tree). Therefore, the
computation required to locate the two nearest clusters is equal to the number of trees multiplied by
the computation required to search for the smallest element of each tree. This yields a total search
time which is O(M log(M \Gamma 1)). After two clusters have been merged a new tree must be created,
and the remaining trees must be updated by deleting the old cluster entries and inserting the
new entry. Together these operations require O(M log M) computation. Therefore, the complete
operation of combining two clusters requires O(M log M) computation, and the total complexity
of the agglomerative clustering procedure is O(M 2 log M ). It should be noted that this algorithm
requires O(M 2 ) memory.
Alternatively, a simpler algorithm which does not assure steepest descent may also be used. This
algorithm, which is shown in Figure 2, starts by randomly ordering the clusters into a linked list. A
cluster is then selected, and the list is searched to locate the nearest neighbor to the selected cluster.
If c(k; l) is negative, the clusters are combined to form a new cluster which replaces the original
two clusters, and the new cluster is selected. Otherwise, a cluster is selected from a member of the
list which has not been previously selected. This algorithm requires only O(M 2 ) computation, and
O(M) memory.
A direct method for clustering the data is to start off with each image block as an individual
cluster. However, the initial number of blocks (ranging from 1024 to 4096 in our experimental
results) is often is too large to be easily handled. In addition, at this stage many cluster pairs result
in a negative value for c(k; l). Therefore, we first separate the data into quantization bins based
on the mean and standard deviation of the sample statistics for a block. This efficiently groups
together blocks which are very similar in their statistics. In all experiments we used 64 quantization
levels for mean and 64 quantization levels for standard deviation for a total of 4096 distinct cluster
types. This coarse grouping reduced the number of clusters to a range between 200 and 500.
After this initial grouping, all merging of clusters is done to strictly reduce the clustering criteria
of (23). Since the initial number of clusters is large, the algorithm of Figure 2 is first used to merge
all clusters with c(k; l) ! 10. From this point on, only the steepest descent algorithm was used.
This strategy is conservative since both these algorithms were found to yield similar and usually
equivalent results when used individually.
An initial assumption of this analysis is that each block contains only a single texture. However,
there will generally be blocks which fall on region boundaries and contain a mixture of textures.
Zhang and Modestino have proposed a method of testing the variance of pixels in an image block
and removing blocks that have excessively high variance [16]. However, in texture segmentation this
method is unreliable since normal image textures are likely to have large variance. The approach
we take is to perform ICM segmentation at the coarsest level of resolution using the estimated
texture parameters. All blocks which fall on a texture boundary are then identified. A boundary
block is defined as any block which differs in texture from one of its four nearest neighbors. All of
these boundary blocks are then removed from their associated clusters. Cluster statistics and ML
parameters are then recalculated and the clusters are merged using the steeped descent algorithm.
The complete algorithm for unsupervised parameter estimation is given below.
Extract texture statistics for each block of the image by using (18).
Group these statistics coarsely by mean and prediction variance in order to generate a tractable
number of clusters.
3 Calculate ML parameter estimates for each cluster.
4 Combine clusters until no more negative values of c(k; l) occur.
5 Calculate the cluster partitions, fA k g M
k=1 , which minimize (23). If any changes have occurred
return to step 3. Otherwise continue.
6 Perform ICM segmentation at the coarsest resolution, and remove from consideration any
blocks which fall on boundaries. Recalculate the ML parameters, and then combine the
remaining clusters until no more negative values of c(k; l) occur.
To verify the performance of the MRS algorithm, computer simulations were done under a
variety of conditions which stressed the performance of the MRS, ICM and SA algorithms. The
computational requirements of the segmentation algorithms were measured using two methods. The
first method counts the number of visits per pixel at full resolution. In the case of MRS, this will
usually be a fractional number since a visit at coarse resolution only represents a fraction of a visit
at full resolution. However, this measure is very conservative for both deterministic algorithms since
most local operations only require the testing of a flag (described at the end of section 3.1). The
number of actual replacement computations per pixels is a more relevant measure of computation
for implementations with significantly fewer computing elements than pixels.
All plots are in terms of the parameter - which represents the cost per unit length of boundaries
between differing regions.
The constants relating - and the parameters fi are chosen to meet two constraints when used in
conjunction with a square two dimensional lattice, S, with unit spacing between samples on both
axes. First, the energy function should be the same for diagonal, horizontal and vertical
boundaries of equal length. Second, for these boundaries the energy function should equal the unit
boundary length multiplied by -. Both these constraints assume that boundaries are long (ignore
end effects), and neither holds for boundaries which fall at angles other than diagonal, horizontal
and vertical.
In all experiments, SA used an annealing schedule of the form
where T (n) is the temperature used in the n th pass. This annealing schedule is uniquely determined
by the starting and ending temperatures and the number of iterations. In all experiments, the SA
algorithm was started with the ML estimate of X as the initial condition.
There is some experimental evidence to believe that varying - as a function of the number of
iterations can improve performance in some applications of SA and ICM. For example, Besag [4]
found that progressively increasing - improved the performance of ICM in certain cases. Also, a
number of studies have found that varying the parameters of the Gibbs distribution during ICM
or SA can improve results when modeling continuously valued fields with discontinuities [10, 12].
However, in the case of continuous valued fields, the parameters were progressively reduced from
initially large values. In any case, the schedule needed for changing these parameters introduces
additional free parameters which generally are image dependent. Therefore, our approach is to fix
the parameter - when using ICM, or SA; and, as described earlier, to fix the parameter - at all
levels of resolution when using MRS.
A series of tests were performed on two simple 64x64 toroidal test patterns shown in Figure
3 to determine the relative performance of ICM, MRS and SA. In all these tests, three levels of
resolution were used, and 100 interactions of SA were performed with 1=T

Figures

4 and 5 plot the mean number of misclassified pixels and the value of the
energy (as defined in (7) ) as a function of - for MRS, ICM and SA. A set of plots is shown for two
different cases each of which uses a test pattern shown in Figure 3 with added Gaussian noise. The
signal-to-noise ratio was defined as
The plots of mean error versus - indicate that the performance of the MRS and SA algorithms
are relatively insensitive to the choice of -. The dynamic range of the energy function plotted in

Figure

5 makes relative differences for pattern 1 difficult to distingish. The curves for SA and MR
fall very close and effectively overlap in the plot.
Three additional tests were done for fixed - to compare the performance of MRS, ICM and SA.
The empirical density function of the error is shown in Figures 6, 7 and 8 together with the mean
percentage error, the mean energy of the segmentation and the mean number of visits per pixel for
each of the three methods. In all these examples, the MRS algorithm requires less computation
than either ICM or SA. The empirical error distributions and the mean energy values indicate that
the performance of MRS is comparable to SA and substantially better than ICM in some cases. In
general, the advantages of MRS are greatest when the signal-to-noise ratio is low. In this case, the
information needed to accurately classify pixels is spread over larger regions.
Segmentation was also performed for six different tests shown in Figures 9 to 14. In all but Figure
unsupervised parameter estimation was used. All test images have a resolution of 256x256,
and all texture models used a three point prediction model using the three
closest pixels in the prediction. In each case, a free boundary was assumed with the parameter value
unsupervised parameter estimation was performed, it was done at the coarsest resolution
used in segmentation. Each example contains segmentations resulting from MRS, ICM and
SA for 100 and 500 iterations, a starting temperature of T and an ending temperature
of T possible the correct segmentation was also presented. Each distinct
texture is displayed as a unique gray level outlined in white.

Tables

summerize the results of the unsupervised estimation algorithm. Table 1 lists
the number of distinct texture types estimated for each of the five experiments shown in Figures 10
to 14. The true number of textures is also listed for the three experiments in which the number of
textures was known. In each case, the estimated number of textures was correct. Table 2 shows
the actual and estimated parameter values for the synthetic image used in Figures 9 and 10. The
actual parameters were those used in generating the synthetic textures.

Figures

9 to 12 used composite images formed by piecing together both natural and synthetic
textures. In these composite images, adjacent textures were smoothly overlapped over a two pixel
boundary in order to avoid a visible cue. In all four tests, four levels of resolution were used
It should be noted that the unsupervised parameter algorithm requires at least one region of each
texture large enough to guarantee that the region has an interior and therefore will be detected.
The pattern for the texture regions has been chosen so this is true with four levels of resolution.

Figures

both used the same synthetic image. This composite image was formed
from synthetic textures generated using the assumed texture model. Figure 9 shows the result
of segmenting the synthetic image using the actual texture parameters, and Figure 10 shows the
result of segmenting the same image using the parameters obtained from the unsupervised parameter
estimation algorithm. Figures 11 and 12 used composite images formed from natural texture types.
In order of their area, the first image contains woolen cloth, raffia and grass; and the second image
contains woolen cloth, grass and beach sand.

Figures

13 and 14 each show the segmentation of an aerial photograph. Three levels of resolution
are used so that smaller regions with distinct textures could be identified. For these images,
the quality of the segmentation depends on the application, and it is possible to vary the amount
of detail in the segmentation by varying the boundary cost -. Nevertheless, a superior algorithm
should more reproducibly classify similar regions to the same texture.

Table

4 lists the value of the energy function resulting from each of the three minimization
algorithms MRS, ICM and SA. The energy of the segmentation is useful since it provides an objective
measure of performance. The computation required for segmentation using MRS, ICM and SA is
shown in Table 3 both in terms of visits per pixel and replacements per pixel. In every case, the
MRS algorithm required less computation than ICM and much less than SA. It should be noted
that the number of replacements per pixel required by the MRS algorithm is relatively independent
of the image segmented.
6 Conclusion
The MRS algorithm proposed in this paper performed better and required less computation than
ICM in all cases tested. Also, the MRS algorithm generally performed comparably to SA but with
substantially less computation. In fact, the computation required for SA ranged between 10 to 1000
times that of MRS depending on the problem and the method of comparison. However, the number
of replacements per pixel was found to be a good basis for comparison, and using this measure SA
typically required two orders of magnitude more computation for a range of examples. The MRS
algorithm seems to yield the most substantial improvement in images when the information per
pixel is low, so that information over larger regions must be combined to correctly segment the
image.
An algorithm for performing unsupervised parameter estimation was proposed for use in conjunction
with the segmentation algorithm. This algorithm is based on a Gaussian AR texture
model and estimates both the number of textures with statistically significant differences and the
parameters of those textures.
The algorithm is composed of three basic operations. The partitioning of data into distinct
textures or clusters, the estimation of the cluster parameters and the agglomeration of similar
clusters to reduce the number of distinct textures. All of these operations are performed in the
context of minimizing a proposed statistical criteria, so the algorithm is guaranteed to converge.
Superfluous textures formed at texture boundaries are removed in a final step.
A

Appendix


In this section we will derive the joint distribution of Y and fA k
k=1 under the assumptions
that the relative order of blocks in each set A k is ignored, and that the number of blocks of each
texture, W k , have a multinomial distribution with unknown parameters, ae k . Using the function
L(\Delta; \Delta) defined in (17), we can rewrite the conditional distribution of Y given the partitioning sets
k=1 explicitly as a function of the statistics, - A k , and the parameters ' k . This yields the
expression
Y
where
One method of conceptualizing the loss of ordering information is to assume that the values in
X (L\Gamma1) and the corresponding statistics, - s , have been reordered with a random mapping that is
uniformly distributed over all possible reorderings. In this case, any unique mapping occurs with
equal probability 1=W ! . If \Pi s is the random reordering, then
If we condition on the knowledge of fW k
k=1 and the value of X (L\Gamma1) , then the probability of any
particular group of partitioning sets is given by the sum of the probabilities of all mappings which
would result in that partitioning. Since there are exactly Q M
unique mappings which result
in any group of partitioning sets, we have that
We next assume that fW k
k=1 have a multinomial distribution with unknown parameters
Y
By combining (28) with (29) the final expression of the joint distribution of Y and fA k
results.
Y
The distribution for Y can now be computed by summing over all possible partitioning sets, or
equivalently, by summing over all possible values of x. To do this, we let the sequence fs
be some ordering of all the points in S. Then we have
xs N =1
xs N =1
Y
ae xs i
expf\Gammal s i (yjx s i )g
Y
ae xs i
expf\Gammal
xs N =1
ae xs N
expf\Gammal s N (yjx s N )g5
Y
ae xs i
expf\Gammal s i (yjx s i )g
where x s is used in place of x (L\Gamma1)
s . By exchanging products and summations, we then obtain
Y
which is now recognizable as the distribution of W independent and identically distributed random
variables with a marginal distribution formed by a mixture of M distributions.



--R

"Modeling and Segmentation of Noisy and Textured Images Using Gibbs Random Fields,"
"Segmentation of Noisy Textured Images Using Simulated Annealing,"
"A Parallel Image Segmentation Algorithm Using Relaxation with Varying Neighborhoods and Its Mapping to Array Processors,"
"On the Statistical Analysis of Dirty Pictures,"
"Spatial Interaction and the Statistical Analysis of Lattice Systems"
"Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images,"
"Advances in Mathematical Models for Image Processing"
"Estimation and Choice of Neighbors in Spatial-Interaction Models of Images,"
"Equations of State Calculations by Fast Computing Machines,"
"Computing Motion Using Analog and Binary Resistive Networks,"
"Segmentation of Textured Images Using a Multiple Resolution Approach,"
"A Multiple Resolution Approach to Regularization,"
"Image Analysis Using Multigrid Relaxation Methods,"
"A New Look at the Statistical Model
"A Model-Fitting Approach to Cluster Validation with Application to Stochastic Model-Based Image Segmentation,"
"Unsupervised Image Segmentation Using A Gaussian Model,"
Pattern Classification and Scene Analysis
Chapman and Hall
"Maximum Likelihood from Incomplete Data via the EM Algorithm,"
"Mixture Densities, Maximum Likelihood and the EM Algorithm,"
"A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains,"
The Art of Computer Programming: Volume
--TR
Image analysis using multigrid relaxation methods
Modeling and segmentation of noisy and textured images using Gibbs random fields
A parallel image segmentation algorithm using relaxation with varying neighborhoods and its mapping
Computing Motion Using Analog and Binary Resistive Networks
The art of computer programming, volume 3

--CTR
Feng Li , Jiaxiong Peng , Xiaojun Zheng, Object-based and semantic image segmentation using MRF, EURASIP Journal on Applied Signal Processing, v.2004 n.1, p.833-840, 1 January 2004
Walter G. Kropatsch , Yll Haxhimusa , Zygmunt Pizlo , George Langs, Vision pyramids that do not grow too high, Pattern Recognition Letters, v.26 n.3, p.319-337, February 2005
Feng Li , Jiaxiong Peng, Double random field models for remote sensing image segmentation, Pattern Recognition Letters, v.25 n.1, p.129-139, 5 January 2004
Zhuowen Tu , Song-Chun Zhu, Image Segmentation by Data-Driven Markov Chain Monte Carlo, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.24 n.5, p.657-673, May 2002
Olivier Alata , Clarisse Ramananjarasoa, Unsupervised textured image segmentation using 2-D quarter plane autoregressive model with four prediction supports, Pattern Recognition Letters, v.26 n.8, p.1069-1081, June 2005
Wei Fu , Guang-Zhong Xing, Edge-oriented spatial interpolation for error concealment of consecutive blocks, Journal of Computer Science and Technology, v.22 n.3, p.494-497, May 2007
Roland Wilson , Chang-Tsun Li, A Class of Discrete Multiresolution Random Fields and Its Application to Image Segmentation, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.25 n.1, p.42-56, January
A. Lorette , X. Descombes , J. Zerubia, Texture Analysis through a Markovian Modelling and FuzzyClassification: Application to Urban Area Extraction fromSatellite Images, International Journal of Computer Vision, v.36 n.3, p.221-236, Feb.-March 2000
Majid Mirmehdi , Maria Petrou, Segmentation of Color Textures, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.22 n.2, p.142-159, February 2000
Nikos Paragios , Rachid Deriche, Geodesic Active Regions and Level Set Methods for Supervised Texture Segmentation, International Journal of Computer Vision, v.46 n.3, p.223-247, February-March 2002
Djamal Boukerroui , Atilla Baskurt , J. Alison Noble , Olivier Basset, Segmentation of ultrasound images: multiresolution 2D and 3D algorithm based on global and local statistics, Pattern Recognition Letters, v.24 n.4-5, p.779-790, February
Alexandra Pizurica , Wilfried Philips , Ignace Lemahieu , Marc Acheroy, The application of Markov random field models to wavelet-based image denoising, Imaging and vision systems: theory, assessment and applications, Nova Science Publishers, Inc., Commack, NY, 2001
Olivier Alata , Christian Olivier, Choice of a 2-D causal autoregressive texture model using information criteria, Pattern Recognition Letters, v.24 n.9-10, p.1191-1201, 01 June
Yue Wang , Tlay Adali , Chi-Ming Lau , Sun-Yuan Kung, Quantitative Analysis of MR Brain Image Sequences by Adaptive Self-Organizing Finite Mixtures, Journal of VLSI Signal Processing Systems, v.18 n.3, p.219-239, April 1, 1998
Murilo P. Almeida, Anisotropic Textures with Arbitrary Orientation, Journal of Mathematical Imaging and Vision, v.7 n.3, p.241-251, June 1997
Ayman El-Baz , Aly A. Farag , Georgy Gimel'farb, Iterative approximation of empirical grey-level distributions for precise segmentation of multimodal images, EURASIP Journal on Applied Signal Processing, v.2005 n.1, p.1969-1983, 1 January 2005
C. C. Reyes-Aldasoro , A. Bhalerao, The Bhattacharyya space for feature selection and its application to texture segmentation, Pattern Recognition, v.39 n.5, p.812-826, May, 2006
statistical behavior;classification;variance;picture processing;posteriori probability;deterministic greedy algorithm;markov random field;statistics;textured images;pattern recognition;probability;parameter estimation;unsupervised texture segmentation;spatial correlation;multiple resolution segmentation;causal Gaussian autoregressive model;coarse resolution;AIC information criterion