A Robust Competitive Clustering Algorithm With Applications in Computer Vision.
AbstractThis paper addresses three major issues associated with conventional partitional clustering, namely, sensitivity to initialization, difficulty in determining the number of clusters, and sensitivity to noise and outliers. The proposed Robust Competitive Agglomeration (RCA) algorithm starts with a large number of clusters to reduce the sensitivity to initialization, and determines the actual number of clusters by a process of competitive agglomeration. Noise immunity is achieved by incorporating concepts from robust statistics into the algorithm. RCA assigns two different sets of weights for each data point: the first set of constrained weights represents degrees of sharing, and is used to create a competitive environment and to generate a fuzzy partition of the data set. The second set corresponds to robust weights, and is used to obtain robust estimates of the cluster prototypes. By choosing an appropriate distance measure in the objective function, RCA can be used to find an unknown number of clusters of various shapes in noisy data sets, as well as to fit an unknown number of parametric models simultaneously. Several examples, such as clustering/mixture decomposition, line/plane fitting, segmentation of range images, and estimation of motion parameters of multiple objects, are shown.
Introduction
Traditional clustering algorithms can be classified into two main categories [1]: hierarchical
and partitional. In hierarchical clustering, the number of clusters need not be specified a
priori, and problems due to initialization and local minima do not arise. However, since
hierarchical methods consider only local neighbors in each step, they cannot incorporate
a priori knowledge about the global shape or size of clusters. As a result, they cannot
always separate overlapping clusters. Moreover, hierarchical clustering is static, and points
committed to a given cluster in the early stages cannot move to a different cluster.
Prototype-based partitional clustering algorithms can be divided into two classes: crisp
(or hard) clustering where each data point belongs to only one cluster, and fuzzy clustering
where every data point belongs to every cluster to a certain degree. Fuzzy clustering algorithms
can deal with overlapping cluster boundaries. Partitional algorithms are dynamic,
and points can move from one cluster to another. They can incorporate knowledge about
the shape or size of clusters by using appropriate prototypes and distance measures. These
algorithms have been extended to detect lines, planes, circles, ellipses, curves and surfaces
[2, 3, 4, 5]. Most partitional approaches use the alternating optimization technique, whose
iterative nature makes them sensitive to initialization and susceptible to local minima. Two
other major drawbacks of the partitional approach are the difficulty in determining the
number of clusters, and the sensitivity to noise and outliers.
In this paper, we describe a new approach called Robust Competitive Agglomeration
(RCA), which combines the advantages of hierarchical and partitional clustering techniques
[6]. RCA determines the "optimum" number of clusters via a process of competitive agglomeration
while knowledge about the global shape of clusters is incorporated via the
use of prototypes. To overcome the sensitivity to outliers, we incorporate concepts from
robust statistics. Overlapping clusters are handled by the use of fuzzy memberships. The
algorithm starts by partitioning the data set into a large number of small clusters which
reduces its sensitivity to initialization. As the algorithm progresses, adjacent clusters compete
for points, and clusters that lose the competition gradually vanish. However, unlike in
traditional hierarchical clustering, points can move from one cluster to another. RCA uses
two different sets of weights (or memberships) for each data point: the first one is a set of
probabilistically constrained memberships that represent degrees of sharing among the clus-
ters. The constraint generates a good partition and introduces competition among clusters.
The second set of memberships is unconstrained or possibilistic [8, 9, 10], and represents
degrees of "typicality" of the points with respect to the clusters. These memberships are
used to obtain robust estimates of the cluster prototypes.
The organization of the rest of the paper is as follows. In section 2, we briefly review
other related approaches. In section 3, we present the RCA algorithm. In section 4, we
illustrate the power and flexibility of RCA to incorporate various distance measures. In
section 5, we describe the application of RCA to segmentation of range images. In section 6,
we formulate a multiple model general linear regression algorithm based on RCA and apply
it to simultaneous estimation of motion parameters of multiple objects. Finally, section 7
contains the conclusions.
Related Work
Most prototype-based partitional clustering algorithms such as K-Means and Fuzzy C-Means
assume that the number of clusters, C, is known. Moreoever, since they use a
least squares criterion, they break down easily (i. e., the prototype parameter estimates
can be arbitrarily wrong [11]) in the presence of noise. The goal of clustering is to identify
clusters in the data set. This implicitly assumes that we have a definition for a valid cluster.
Thus, the idea of break down [11] can be extended to the clustering domain via the use
of validity [12]. When the number of clusters, C, is known, the ideal cluster breaks down
only when the outliers form a valid cluster with a cardinality higher than the cardinality,
min , of the smallest good cluster. This gives us the theoretical breakdown point of N min =N ,
where N is the number of points in the data set. Recent solutions to robust clustering when
C is known can be divided into two categories. In the first category are algorithms that
are derived by modifying the objective function of FCM [13, 14, 10]. These algorithms are
still sensitive to initialization and other parameters [12]. The algorithms in second category
incorporate techniques from robust statistics explicitly into their objective functions. A
notable non-fuzzy clustering algorithms in this category is the K-Medoids algorithm [15].
Bobrowski and Bezdek [16] proposed an L 1 -norm-based fuzzy clustering algorithm which
also falls into this category. However, there is no mention of robustness in this paper. A
variation of this algorithm that is motivated by robustness can be found in [17]. Another
early fuzzy clustering algorithm (on which RCA is based) is the Robust C-Prototypes (RCP)
algorithm [18], which uses the M-estimator [19]. The Fuzzy Trimmed C Prototypes
algorithm [20] uses the least trimmed squares estimator [21], the Robust Fuzzy C Means
(RFCM) algorithm [22] again uses the M-estimator in a different way, and the Fuzzy C
Least Median of Squares (FCLMS) algorithm [23] uses the least median of squares estimator
[21]. FTCP and FCLMS can achieve the theoretical breakdown point of N min =N with a
trivial modification to their objective functions. However, in theory, they both require an
exhaustive search. To reduce the computational complexty, a heuristic search is used in
[20]and a genetic search is used in [23].
When C is unknown, one way to state the clustering problem is: find all the valid clusters
in the data set (see [12] for a more precise definition). In this case, the ideal algorithm will
not break down because it will identify all the "good" clusters correctly (say by exhaustive
search), in addtion to some spurious ones. An alternative way to state the problem is:
identify only all the valid clusters formed by the good data. In this case, the ideal algorithm
will break down when the outliers form a valid cluster, giving us the breakdown point of
N minval =N , where N minval is the minimum number of points required to form a valid cluster.
Note that a given clustering algorithm may not achieve these theoretical breakdown points.
The traditional approach to determining C is to evaluate a certain global validity measure
of the C-partition for a range of C values, and then pick the value of C that optimizes
the validity measure [25, 1, 26, 27]. An alternative is to perform progressive clustering
[28, 27, 5], where clustering is initially performed with an overspecified number of clusters.
After convergence, spurious clusters are eliminated, compatible clusters are merged, and
"good" clusters are identified. Another variation of progressive clustering extracts one cluster
at a time [29, 30]. These approaches are either computationally expensive, or rely on validity
measures (global or individual) which can be difficult to devise. Robust approaches to
clustering when C is unknown treat the data as a mixture of components, and use a robust
estimator to estimate the parameters of each component. The Generalized MVE (GMVE)
[29] which is based on the Minimum Volume Ellipsoid estimator [21], the Model Fitting (MF)
algorithm [31], and the Possibilistic Gaussian Mixture Decomposition (PGMD) algorithm
[30] are some examples. In the above approaches, the data set is classified into a set of
"inliers", i.e., points belonging to a cluster, and a set of "outliers". Since the set of outliers
includes points from other clusters, the proportion of outliers can be very high. Therefore,
even the use of a robust estimaor with the theoretical-best breakdown point of 50% is not
sufficient to make these algorithms highly robust. To overcome this problem, these algorithms
consider the "validity" of the cluster formed by the inliers, and try to extract every valid
cluster in the data set. In order to guarantee a good solution, the GMVE and PGMD use
many random initializations. Cooperative Robust Estimation (CRE) [32] and MINPRAN
[33] are two other robust model-fitting approaches that fall into this category. The CRE
algorithm attempts to overcome the low breakdown point of M-estimators by initializing a
large number of hypotheses and then selecting a subset of the initial hypotheses based on the
Minimum Description Length (MDL) criterion. The CRE technique assumes that the scale
(' in [32]) is known. MINPRAN assumes that the outliers are randomly distributed within
the dynamic range of the sensor, and the noise (outlier) distribution is known. Because of
these assumptions, CRE and MINPRAN do not easily extend to the clustering domain. If
the data is expected to have multiple curves, MINPRAN seeks one curve/surface at a time.
In [12] the relation between the above progressive approaches and other robust clustering
algorithms are explored.
When the clusters overlap, the idea of extracting them in a serial fashion will not work.
Removing one cluster may partially destroy the structure of other clusters, or we might get
"bridging fits" [33]. Fig. 2(a) shows one such noisy data set with two crossing clusters.
The algorithm we propose is designed to overcome this drawback. Moreover, all the current
algorithms use hard finite rejection [34], i.e., points within an inlier bound are given a
weight of 1, and points outside the bound are given a weight of zero. This means that these
algorithms do not handle the "region of doubt" [21] very well. To overcome this problem,
we use smooth [34, 21] or fuzzy rejection, where the weight function drops to zero gradually.
3 The Robust Competitive Agglomeration (RCA) algorith

3.1 Algorithm Development
be a set of N vectors in an n-dimensional feature space with
coordinate axis labels C-tuple of prototypes
each of which characterizes one of the C clusters. Each fi i consists of a set of parameters.
The Fuzzy C-Means algorithm [2] minimizes:
subject to
In (1), d 2
ij represents the distance of feature vector x j from prototype fi i , represents the
degree to which x j belongs to cluster i, is a C \ThetaN matrix called the constrained
fuzzy C-partition matrix, and m 2 [0; 1) is known as the fuzzifier. Jm , which is essentially
the sum of (fuzzy) intra-cluster distances, has a monotonic tendency, and has the minimum
value of zero when C=N . Therefore, it is not useful for the automatic determination of C.
To overcome this drawback, we add a second regularization term to prevent overfitting the
data set with too many prototypes. The resulting objective function JA is:
which is minimized subject to the constraint in (2). In (3), the second term is the negative
of the sum of the squares of the cardinalities of the clusters, and is minimized when the
cardinality of one of the clusters is N and the rest of the clusters are empty. With a proper
choice of ff, we can balance the two terms to find a solution for C. JA is still not robust,
since the first term is a Least Squares objective function. Therefore, we robustify JA to yield
the objective function for the proposed RCA algorithm as follows:
In (4), ae i () is a robust loss function associated with cluster i, and w
represents the "typicality" of point x j with respect to cluster i. The function ae i () corresponds
to the loss function used in M-estimators of robust statistics and w i () represents the weight
function of an equivalent W-estimator (see, [11], for example). This particular choice for
robustification is motivated by the need to keep the computational complexity low. The loss
function reduces the effect of outliers on the first term, and the weight function discounts
outliers while computing the cardinalities. By selecting d ij and the ff prudently, JR can be
used to find compact clusters of various types while partitioning the data set into a minimal
number of clusters.
To minimize JR with respect to the prototype parameters, we fix U and set the derivative
of JR with respect to fi i to zero, i.e.,
0: (5)
Further simplification of (5) depends on ae i () and d ij . Since the distance measure is application
dependent, we will return to this issue in Section 4. To minimize (4) with respect to U
subject to (2), we apply Lagrange multipliers and obtain
We then fix B and solve
st
st ae s (d 2
st
Equations (7) and (2) represent a set of of N \ThetaC +N linear equations with N \ThetaC +N unknowns
st , and - t ). A computationally simple solution can be obtained by computing the
using the memberships from the previous iteration. This yields:
st
2ff \Theta (
st
Solving for - t using (8) and (2), and substituting in (8), we obtain the following update
equation for the membership u st of feature point x t in cluster
st
st
ff
ae s (d 2
st
st
st
where u RR
st is the degree to which cluster s shares x t (computed using robust distances), and
st is a signed bias term which depends on the difference between the robust cardinality,
of the cluster of interest and the weighted average of cardinalities
kt
. C
kt
The bias term, u Bias
st , is positive(negative) for clusters with cardinality higher(lower) than
average, and hence the membership of x t in such clusters will appreciate(depreciate). When
a feature point x j is close to only one cluster (say cluster i), and far from other clusters,
we have N i - N j , or u Bias
implying no competition. On the other hand, if a point
is roughly equidistant from several clusters, these clusters will compete for this point based
on cardinality. When the cardinality of a cluster drops below a threshold, we discard the
cluster, and update the number of clusters.
It is possible for u ij to become negative if N i is very small and point x j is close to other
dense clusters. In this case, it is safe to set u ij to zero. It is also possible for u ij to become
larger than 1 if N i is very large and feature point x j is close to other low cardinality clusters.
In this case it is clipped to 1. This practice is customary in optimization theory.
The process of agglomeration, controlled by ff, should be slow in the beginning to encourage
the formation of small clusters. Then it should be increased gradually to promote
agglomeration. After a few iterations, when the number of clusters becomes close to the
"optimum", the value of ff should again decay slowly to allow the algorithm to converge.
Therefore an appropriate choice of ff in iteration k is.
ij
In (10), ff and j are functions of the iteration number k, and the superscript used
on
ij , and w ij to denote their values in iteration k \Gamma 1. A good choice for j is
is the initial value, - is the time constant, and k 0 is the iteration number at which
j starts to decrease. In all examples presented in this paper (except in section 5 where
these parameters were fine-tuned for best performance), we choose
proper initialization, these values are reasonable regardless of the application.
Initialization issues are discussed in section 7.
3.2 Choice of the weight function
In curve/surface fitting or linear regression, it is reasonable to assume that the residuals
have a symmetric distribution about zero. Therefore, we choose Tukey's biweight function
[11] given by
if jr
where r
ij stands for the normalized residual defined as:
r
c \Theta MAD i
In (12)-(14), r ij is the residual of the j th point with respect to the i th cluster, Med i is the
median of the residuals of the i th cluster, and MAD is the median of absolute deviations [11]
of the i th cluster. In other words, in each iteration, the data set X is crisply partitioned into
Med i and MAD i are estimated for each cluster.
When distances (rather than residuals) are used, the symmetric distribution assumption
does not hold. We suggest a monotonically non-increasing weight function w
[0; 1] such that w i (d 2 is a constant, and T i and S i are given
by
Choosing w results in the following weight function:
The corresponding loss function can be shown to be
In (17) K i is an integration constant used to make all ae i () reach the same maximum value.
This choice ensures that all noise points will have the same membership value in all clusters.
Fig 1 shows the plot of the weight function and the corresponding loss function.
In (14), (16), and (17), c is a tuning constant [11] which is normally chosen to be between
4 and 12. When c is large, many outliers will have small nonzero weights, thus affecting the
parameter estimates. On the other hand, if c is small, only a subset of the data points will
be visible to the estimation process, making convergence to a local minimum more likely. As
a compromise, we start the estimation process with a large value of c, and then decrease it
gradually as function of the iteration number (k), i.e.,
with c 0 =12, c min =4, and \Deltac=1.
The RCA algorithm is summarized below.
Fix the maximum number of clusters
Initialize the prototype parameters, and set
Repeat
Compute
Estimate T i and S i by using (15);
Update the weights w ij by using (13) or (16);
Update ff(k) by using (10);
Update the partition matrix U (k) by using (9);
Compute the robust cardinality N
If
Update the number of clusters C;
Update the tuning factor c by using (18);
Update the prototype parameters;
Until
prototype parameters stabilize
4 Examples of Distance Measures
As mentioned in section 3.1, RCA can be used with a variety of distance measures depending
on the nature of the application. In this section, we discuss distance measures suitable for
ellipsoidal clusters and hyperplanes.
4.1 Detection of Ellipsoidal Clusters
To detect ellipsoidal clusters in a data set, we use the following distance measure [35, 36].
In (19), c i is the center of cluster fi i , and C i is its covariance matrix. (See [37] for an
interpretation of d 2
Cij .) Using (5), it can be shown that the update equations for the centers
c i and the covariance matrices C i are
If we assume C reduces to the Euclidean distance. This simplified version
can be used when the clusters are expected to be spherical.
Fig. 3 illustrates RCA using d 2
Cij . Fig. 3(a) shows a synthetic Gaussian mixture consisting
of 4 clusters of various sizes and orientations. Uniformly distributed noise constituting
40% of the total points was added to the data set. Fig. 3(b) shows the initial 20 prototypes
superimposed on the data set, where "+" signs indicate the cluster centers, and the ellipses
enclose points with a Mahalanobis distance less than 9. These prototypes were obtained by
running the G-K algorithm [36] for 5 iterations. After 2 iterations of RCA, 9 empty clusters
are discarded (see Fig. 3(c)). The number of clusters is reduced to 6 after 3 iterations, and
to 4 after 4 iterations. The final result after a total of 10 iterations is shown in Fig. 3(d).
To illustrate the ability of RCA to handle non-uniform noise, Fig. 4 shows the result of
RCA on a data set containing Gaussian clusters with roughly 25% noise. To illustrate the
ability of the RCA algorithm to detect overlapping clusters, in Fig. 2(b) we show the result
of RCA on the data set in Fig. 2(a). The algorithm converged in 10 iterations.
4.2 Detection of Linear Clusters
To detect clusters that resemble lines or planes, we use a generalization of the distance
measure proposed in [3, 2]. This distance is given by
where e ik is the k th unit eigenvector of the covariance matrix C i . The eigenvectors are
assumed to be arranged in ascending order of the corresponding eigenvalues. The value of
- ik in (22) is chosen dynamically in every iteration to be - is the k th
eigenvalue of C i . It can be shown that for the distance measure in (22), the update equations
for c i and C i are given by (20) and (21) respectively.
Fig. 5(a) shows an image consisting of 10 line segments in a noisy background. Fig. 5(b)
shows the 20 initial prototypes obtained by running the AFC algorithm [3] for 5 iterations.
After 2 iterations of RCA, the number of clusters drops to 15 as shown in Fig. 5(c). After
9 iterations, the number of clusters reduces to the "optimal" number and the algorithm
converges after a total of 12 iterations. The final result is shown in Fig. 5(d).
5 Application to Range Image Segmentation
5.1 Planar Range Image Segmentation
Since planar surface patches can be modeled by flat ellipsoids, the distance measure d 2
in (19) can also be used to find the optimal number of planar patches. To avoid missing
tiny surfaces, we start by dividing the image into non-overlapping windows of sizes W s \ThetaW s .
Then, we apply RCA in each window with to estimate the optimal number of
planar patches within the window. Finally, we pool the resulting (say M) prototypes to
initialize the RCA algorithm with C=M . Because of the nature of d 2
Cij , planar surfaces
with non-convex shapes may be approximated by several planar patches, or several spatially
disconnected planar patches may be approximated by a single cluster. Therefore, after RCA
converges, we merge compatible clusters [27] that are adjacent. We then perform connected
component labeling on each cluster, and assign different labels to disjoint regions.
The above RCA-based algorithm was tested on two standard data sets, ABW data set
and perceptron data set, that were created for bench-marking range image segmentation
algorithms [38]. Each set contains 40 images of size 512\Theta512, and has been randomly divided
into a 10-image training set and a 30-image testing set. We use the performance measures
developed by Hoover et al. [38] to evaluate the performance of RCA. These measures rely
on comparing the Machine Segmented (MS) image and the Ground Truth (GT) image, and
classify the regions into one of the 5 categories: correct detection, over-segmentation, under-
segmentation, missed, and noise. The accuracy of the segmentation is quatified by computing
the average and standard deviation of the differences between the angles made by all pairs
of adjacent regions that are instances of correct detection in the MS and GT images. The
above data sets and performance measures have been used in [38] to compare the University
of South Florida (USF), University of Edinburgh (UE), Washington State University (WSU),
and University of Bern (UB) segmentation algorithms. Here, we will reproduce the same set
of experiments and include the RCA algorithm in the comparison.
In the training phase, we fine-tuned the parameters of RCA as follows: window size used
in the initialization W initial number of prototypes in each window C
(j (see (10)). These parameters are optimal for both ABW and Perceptron
data sets. Since the Perceptron data is more noisy, we use c 4, and for the ABW
data, c 8. Also, to reduce computations, all images were subsampled in the x and y
directions by a factor of 3. These parameters are all then fixed in the testing phase.
Fig 6(a) shows the intensity image of one of the ABW test images. The segmented range
image is shown in Fig 6(b). The shaded gray regions correspond to background points that
are ignored during segmentation. Fig. 7 shows an example from the Perceptron data set.
As in [38], we compute the performance metrics of the five segmentation algorithms while
varying the compare tool tolerance from 51% to 95%. Due to space limitation, we only show
plots of the correct detection measure (Fig 8). The performance measures using an 80%
compare tolerance for all five segmenters are listed in Table 1 for the ABW data and Table
2 for the Perceptron data. RCA compares very well with the best segmenters.
Among the 5 planar surface segmenters in the comparison, UE, WSU, and RCA have the
capability to segment curved surfaces. RCA has the additional advantage that it can handle
irregularly spaced sparse data as well (e.g. range data computed from stereo methods).
5.2 Quadric Range Image Segmentation
Let the i th prototype fi i , represented by the parameter vector p i , define the equation of a
quadric surface as p T
and is a 3-D point. Since the exact distance from a point x j to a quadric surface
fi i has no closed-form expression, we use the approximate distance [39, 40] given by
is the Jacobian of q evaluated at x j . To avoid the all-zero trivial solution for
the following constraint may be chosen [39]
Starting from (5), it can be shown that the use of d Aij leads to a solution of p i based on the
following generalized eigenvector problem: F i
and G
To obtain a reliable initialization, we divide the image into small non-overlapping win-
dows, and apply RCA in each window with C=1. Finally, we pool the resulting prototype
parameters to initialize the RCA algorithm. Initially, there might be several initial prototypes
corresponding to the same surface. However, due to competition, only one of these
surfaces will survive.
The examples used in this section consist of some 240\Theta240 real and some synthetic range
images 1 . A sampling rate of 3 in the x and y directions was used to reduce computations.
were used to estimate the initial prototypes. Fig. 9(a) shows a synthetic
range image of a plastic pipe. Fig. 9(b) shows the initial 36 surface patches. These patches
were generated after assigning each point to the nearest prototype. Fig. 9(c) shows the final
results, where each each surface is displayed with a different gray value, and the boundaries
are shown in black. Fig. 10(a) shows a real range image of three plastic pipes of different
sizes and orientations. The final results of the RCA algorithm consisting of the correctly
identified surfaces are shown in Fig. 10(b).
To test the robustness of RCA, Gaussian noise (with oe=4) was added to the image in
Fig. 9(a), and about 10% of the data points were randomly altered to become outliers. The
These images were obtained from Michigan State University and Washington State University via anonymous
ftp.
results are shown in Fig. 11, where noise points (i.e. points with zero weight (w ij ) in all
clusters) are shown in black.
6 Estimation of Multiple Motion Groups and Segmen-
tation
In this section, we show how RCA can be used to perform multiple model linear regression,
and apply it to estimation of the motion parameters of multiple motion groups.
6.1 General Linear Regression
The General Linear Regression (GLR) [41] for solving a set of homogeneous equations for
motion parameters can be written as: is the design
matrix with x is the parameter vector, and
is the residual vector. Since the system is homogeneous, we can fix fi
and reformulate the GLR model as: N-dimensional
vector with every component equal to 1, can be
solved by the least squares minimization: min fi
with the solution:
least squares is very sensitive to noise. An alternative is
the weighted least squares: min fi
with the solution: fi
If a data set contains multiple models, the GLR model must be applied repetetively to
extract one model at a time. This approach is computationally expensive, requires models to
be well separated, needs a high breakdown estimator (since while extracting the i th model, all
other models are considered as outliers), and is sensitive to initialization. To deal with these
problems, we propose the Multiple-Model General Linear Regression (MMGLR) method,
which allows the simultaneous estimation of an unknown number of models.
6.2 Multiple-Model General Linear Regression
Let the i th model with the parameter vector fi represented by
is the residual corresponding to the j th data vector in the i th model. MMGRL
minimizes (4) (where d 2
ij is replaced by r 2
subject to the constraint in (2). Solving (5)
corresponding to this situation leads to @
iN ). The resulting update equation for the
parameters is:
In linear regression, it is customary to use the studentized residuals r
where h jj is the j th diagonal element of the hat matrix Huang et al.
[41] showed that the corresponding hat matrix for GLR is H  To
extend this principle to the MMGLR, we compute C hat matrices (i. e., one per model), as
The residuals can be normalized as r
jj . However, this normalization introduces
a bias towards noise points (w ij - 0) or points belonging to other models In
this case h  (i)
jj - 0, and hence no normalization takes place. Also, residuals will be inflated
for points which are typical of the i th model since they are divided by a factor smaller than
one. Therefore, we modify the normalization process as follows:
r
otherwise
where h
In other words, points which are known to be atypical of the
th model, are forced to receive the maximum possible inflation factor.
MMGLR can be used to estimate the motion parameters of multiple objects in the same
scene. The instantaneous velocity -
p(t) of a point located on the surface of
a translating object rotating with an instantaneous angular velocity
is characterized by -
is a vector involving
translation. Let (X(t); Y (t)) be the 2-D prespective projection of p(t) onto the
image plane at Z=1, and let (u(t); v(t)) denote its projective instantaneous velocity. Motion
estimation consists of solving for ! and k using a set of N observations
their corresponding . This can be done by solving
N ]; a
Once h has been determined, the motion parameters
and k can be easily obtained [42]. Since h is 9-dimensional and represents a set
of homogeneous equations, we need only 8 observations to solve for the optical flow [42].
When a scene consists of C independently moving objects, the motion of each object can
be characterized by a different vector h i . In this situation, we need to solve Ah
solves this set of equations where X and fi i correspond to A and h i
respectively. It finds C automatically.
MMGLR requires an overspecified number (C) of initial parameter estimates. We obtain
each one of these estimates by solving on a randomly selected subset of 8 observa-
tions. These C estimates are then pooled together to initialize the MMGLR algorithm. To
ensure a reliable result, the initial number of models C needs to be high. However, since C decreases
drastically in the subsequent iterations, this method is still efficient. Since MMGLR
allows points to move from one model to another, and since fuzzy rejection allows points to
change from inliers to outliers and vice versa smoothly, we can afford to use a smaller number
of initializations than algorithms based on hard rejection. In both experiments described in
this subsection, we use C=50.
Fig. 12(a) shows a synthetic 3-D scene consisting of 4 touching rigid objects, each undergoing
a motion with different rotational and translational velocities. Fig. 12(a) displays
the subsampled and scaled true optic flow field. We contaminated this optic flow field with
Gaussian noise (SNR=70), and additionally altered 20% of the observations randomly to
make them outliers. The resulting optic flow field is shown in Fig. 12(b). MMGLR succeeds
in determining the correct number of motion groups in the scene. It also estimates their motion
parameters accurately, as shown in Table 3. Fig. 12(c) shows the segmented optic flow
field where each motion group is represented by a different symbol. The correctly identified
outliers (points having zero weight w ij in all models) are shown as black dots in Fig. 12(c).
The recovered optic flow field is shown in Fig. 12(d).
Figs. 13(a) and 13(b) show two 512\Theta512 subimages of the 13 th and 14 th frames in
a motion sequence [43] containing a moving truck. In this experiment, the background
motion (due to camera panning) is treated as another motion group to create a multiple
motion scenario. We selected target points on the vehicle, and another 30 points from
the background. The matches of these 60 points were computed using a robust matching
algorithm [44] and verified manually. To illustrate the robustness of MMGLR, we added
another 10 target points with erroneous matches. All 70 points are marked '+' in Figs. 13(a)
and 13(b). The target points and their matches were first converted from pixel coordinates
to image coordinates, and then calibrated [43]. Finally, all target points were integrated
to form the mixture data set f(X is the image
coordinates of the i th target point in the 13 th frame, and displacement vector.
The "ground truth" for the vehicle motion is unknown. Also, since the rotation angle
of the truck is too small (about 5 could not be estimated reliably using two-view point
correspondence and three-view line correspondence algorithms [45]. Since we are testing the
robustness of MMGLR and its ability to detect multiple models, and not the performance of
the linear optic flow algorithm, we compare our results with those obtained when the linear
optic flow algorithm is supplied with the correct data subset for each motion (see Table 4).
MMGLR was first run with only the 60 good target points, and then with the added outliers.
In both cases, the algorithm was able to detect the correct number of motion groups (=2)
and estimate their parameters correctly. Fig. 14 shows the partition of the optic flow field
where the two motion groups and the detected outliers are denoted by different symbols.
7 Discussion and Conclusions
7.1 General Comments
RCA is an attempt at addressing the three main issues of partitional clustering algorithms
(the difficulty in determining the number of clusters, sensitivity to initialization, and sensitivity
to outliers), without sacrificing computational efficiency. RCA minimizes a fuzzy
objective function in order to handle overlapping clusters. Constrained fuzzy memberships
are used to create a competitive environment that promotes the growth of "good" clusters.
Possibilistic memberships [10]are used to obtain robust estimates of the prototype parame-
ters. Concepts from robust statistics have been incorporated into RCA to make it insensitive
to outliers. To handle the region of doubt, and to reduce the sensitivity to initialization,
RCA uses soft finite rejection. The agglomerative property makes it relatively insensitive to
initialization and local minima effects. By using suitable distance measures, we can apply
this algorithm to solve many computer vision problems. The choice of ff in (10) is quite
critical to the algorithm. However, ff can be chosen by trial and error to produce stable results
for a given application. The variety of examples presented in this paper show that this
is possible, and that RCA can provide robust estimates of the prototype parameters even
when the clusters vary significantly in size and shape, and the data set is contaminated.
7.2 Computational Complexity
The RCA algorithm has a computational complexity similar to that of FCM [2], which is
O(NC) in each iteration. Here, N is the number of data points, and C is the number of
clusters. However, additional time is required to estimate the weight function w(d 2 ) which
requires us to compute the median of the squared distances twice (first to compute the
median, then to compute the MAD). The median of a data set can be computed iteratively
using
This procedure converges in O(log N) passes through the data set. Since the distribution of
the squared distances does not change significantly in one iteration, this procedure converges
even faster when the median of the previous iteration is used to initialize the computation
of the median of the current iteration. Thus, the overall complexity can be estimated as
O(N log N +NC) per iteration, or O(NK(logN +C)), where K is the number of iterations.
It is to be noted that the value of C varies from C max to C final . Except for the application
to motion analysis, in all other cases we use a standard algorithm such as FCM to initialize
RCA. Therefore, the initialization overhead is O(NkC max ), where k is a small (- 5) integer.
7.3 Breakdown Issues
As discussed in section 2, when C is known, the breakdown point is N min =N , and when C
is unknown, the breakdown is either undefined or N minv al =N . These results were derived
by the use of validity, and an ideal clustering algorithm would use a validity measure and
an expensive exhaustive search to achieve this level of robustness [12]. However, validity
measures are hard to define in practice unless the distribution of the good points is known.
Moreover, deviations from the assumed distribution can occur with widely varying degrees
in real applications, and it is hard to choose thresholds when their optimal values can vary
widely among different data sets, and even among clusters in the same data set.
RCA is a general purpose algorithm that attempts to achieve robustness with reasonable
computational complexity. This is the rationale behind the choice of the M-estimator to
robustify RCA. This choice limits the breakdown point of RCA to 1
, where p is the dimensionality
of the parameter vector to be estimated. However, since RCA starts with a large
mber of initial prototypes, it is possible to increase its robustness under certain conditions.
RCA uses the initial prototypes to generate a partition. The algorithm consists of updating
the weight function for each component of the partition, then updating the memberships,
and then finally updating the prototypes. This process is repeated until convergence. Since
the weight function uses the median and MAD, it can tolerate up to 50% noise points (within
the component) provided it starts with a good initialization.
Let there be C actual clusters. Let the good points from the k th actual cluster be given the
label "k", let the noise points be labeled "0". Let the (hard) partition
corresponding to the C max initial prototypes be labeled as follows. If a given component
has only noise points, it is labeled "0", otherwise it is labeled "i", where i is the label of
the majority of good points in the component. Let P max
i denote the largest component with
the label i. For the RCA algorithm to give robust results, we require an initialization that
satisfies the following conditions. (i) There exists at least one component that has the label i,
for all The prototype corresponding to P max
i is a good point of the i th actual
cluster. (iii) The largest component labeled "0" is smaller than P
contains more than 50% of points labeled "i". Since the cluster region by definition is denser
than the noise region, by using a sufficiently large number of prototypes, it is usually possible
to achieve an initialization to meet these conditions in practice. Initial prototypes placed in
the cluster region will naturally have larger cardinalities and those in the noise region will
have smaller ones. Conditions (i)-(iv) need to be satisfied in the following iterations as well,
to guarantee that the algorithm will converge to a correct result. However, since cardinalities
are replaced by robust cardinalities in the subsequent iterations, it becomes easier to satisfy
these conditions. When the components coalesce and form the final result, each noise point
will be crisply assigned to one of the components while computing the weight function. In the
worst case, all noise points can be assigned to the smallest cluster. Therefore, conditions (iii)
and (iv) above translate to the requirement that the number of noise points be smaller than
the cardinality of the smallest cluster. Thus, when (i)-(iv) are satisfied, RCA can achieve the
theoretical breakdown point. A similar discussion applies to non-point prototypes as well,
with minor modifications. In this case, each initial prototype can be generated with n data
points, where n is the number of parameters in the prototype.
7.4 Initialization Issues
From the above discussion, it is clear that initialization plays a very important role in the
RCA algorithm. The initialization procedure necessarily varies with the type of prototypes
used, the distance measure used, the type of data, and finally the application. We now
outline some guidelines for initialization.
We can compute a theoretical value for the initial number of clusters, C max , as follows.
Let there be C exp number of actual clusters expected in the data set, let N i denote the
cardinality of cluster i, and let n be the number of points required to generate a prototype.
If we randomly pick n points to generate a prototype, then the probability p that we pick
C exp good prototypes, one from each cluster, is given by
QCexp
. If this selection
is repeated K times, the probability that one of these selections generates good prototypes
for all C exp clusters is given by P . For a given value of P g , we can compute
the value of K as,
log (1\Gammap)
e, and C max can be estimated as C This
value of C max grows exponentially with C exp and n, and therefore is unrealistic.
In practice, an existing clustering algorithm (such as FCM [2], GK [36], AFC [3]) can be
used for initialization. At the end of such an initialization, although not all C max prototypes
are expected to be good, we can assume that each of the C exp clusters has a fairly high
probability, P init
i , of being represented by one of the C max initial prototypes. For example,
consider the case of finding lines in a 2-D data set, i.e. 2. If there are N total points,
there are N(N possible ways to pick a pair of points, and hence N(N possible
random initializations for a line. However, most of these initializations involve points that are
far away from each other and constitute poor initializations. On the other hand, an algorithm
such as AFC will use only nearby points, and the probability that two nearby points belong
to the same line is high. If the data set is an image, then by dividing the image into small
windows and applying a conventional clustering algorithm with a suitable number of clusters
in each window can dramatically increase the value of P init . The probability that all C exp
clusters are represented by the initialization is given by
i . In this case, a much
smaller number of initial clusters will suffice.
Based on the above discussion, we suggest the following rules of thumb. For general clus-
choose C max - N
10\Lambdan
, and use a simple clustering algorithm (such as FCM) to generate
the initial prototypes. Since good points are by definition in dense regions, this initialization
can be expected to meet the conditions discussed in the previous sub-section. The case
of plane and surface fitting can be handled by dividing the image into small windows and
applying a suitable clustering algorithm in each window. In the case of regression, the above
initialization techniques are no longer applicable. Hence, we use a random sampling procedure
to generate the prototypes. Because of this randomness, we require a larger value for
C max . In our applications, we set C max - N

Acknowledgments

The authors would like to thank the anonymous reviewers for their valuable comments.
This work was partially supported by a grant from the Office of Naval Research (N00014-
96-1-0439).



--R

Algorithms for Clustering Data
Pattern Recognition with Fuzzy Objective Function Algorithms
"Use of the adaptive fuzzy clustering algorithm to detect lines in digital images,"
"Adaptive fuzzy c-shells clustering and detection of ellipses,"
"Fuzzy and possibilistic shell clustering algorithms and their application to boundary detection and surface approximation,"
"A robust clustering algorithm based on competitive agglomeration and soft rejection of outliers,"
"Clustering by competitive agglomeration,"
"Fuzzy sets as a basis for a theory of possibility,"
Possibility Theory: An Approach to Computerized Processing of Uncertainty
"A possibilistic approach to clustering,"
Robust Statistics the Approach Based on Influence Functions

"Fuzzy clustering and robust estimation,"

Finding Groups in Data: An Introduction to Cluster Analysis
"c-means clustering with the l 1 and l 1 norms,"
"The fuzzy median and the fuzzy mad,"
"A robust clustering algorithm based on the m- estimator,"
Robust Statistics

John Wiley
"Fuzzy and robust formulations of maximum- likelihood-based gaussian mixture decomposition,"
"A genetic algorithm for robust clustering based on a fuzzy least median of squares criterion,"

"Some new indices for cluster validity,"
"Unsupervised optimal fuzzy clustering,"
"Fitting an unknown number of lines and planes to image data through compatible cluster merging,"
"Progressive fuzzy clustering algorithms for characteristic shape recognition,"
"Robust clustering with applications in computer vision,"
"Gaussian mixture density modeling, decomposition and applications,"
"A highly robust estimator through partially likelihood function modeling and its application in computer vision,"
"Cooperative robust estimation using layers of support,"
"Minpran: A new robust estimator for computer vision,"


"Fuzzy clustering with a fuzzy covariance matrix,"
"Fuzzy and possibilistic clustering methods for computer vision,"
"An experimental comparison of range image segmentation algorithms,"
"Estimation of planar curves, surfaces, and nonplanar space curves defined by implicit equations with application to edge and range image segmentation,"
"A comparison of fuzzy shell-clustering methods for the detection of ellipses,"
"Optic flow field segmentation and motion estimation using a robust genetic partitioning algorithm,"
"A simplified linear optic flow-motion algorithm,"
"A sequence of stereo image data of a moving vehicle in an outdoor scene,"
"A robust technique for matching two uncalibrated images through the recovery of the unknown epipolar geometry,"
"Vehicle-type motion estimation from multi-frame images,"
--TR

--CTR
Ujjwal Maulik , Sanghamitra Bandyopadhyay, Performance Evaluation of Some Clustering Algorithms and Validity Indices, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.24 n.12, p.1650-1654, December 2002
Raffaele Cappelli , Dario Maio , Davide Maltoni, Multispace KL for Pattern Representation and Classification, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.23 n.9, p.977-996, September 2001
Sanghamitra Bandyopadhyay , Ujjwal Maulik, An evolutionary technique based on K-means algorithm for optimal clustering in RN, Information SciencesApplications: An International Journal, v.146 n.1-4, p.221-237, October 2002
Ana L. N. Fred , Anil K. Jain, Combining Multiple Clusterings Using Evidence Accumulation, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.27 n.6, p.835-850, June 2005
M. Scionti , J. P. Lanslots, Stabilisation diagrams: pole identification using fuzzy clustering techniques, Advances in Engineering Software, v.36 n.11-12, p.768-779, November 2005
Miin-Shen Yang , Kuo-Lung Wu, A Similarity-Based Robust Clustering Method, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.26 n.4, p.434-448, April 2004
Cheng-Fa Tsai , Chun-Wei Tsai , Han-Chang Wu , Tzer Yang, ACODF: a novel data clustering approach for data mining in large databases, Journal of Systems and Software, v.73 n.1, p.133-145, September 2004
Bogdan Gabrys, Agglomerative Learning Algorithms for General Fuzzy Min-Max Neural Network, Journal of VLSI Signal Processing Systems, v.32 n.1-2, p.67-82, August-September 2002
A. K. Qin , P. N. Suganthan, Robust growing neural gas algorithm with application in cluster analysis, Neural Networks, v.17 n.8-9, p.1135-1148, October/November 2004
Ana L. N. Fred , Jos M. N. Leito, A New Cluster Isolation Criterion Based on Dissimilarity Increments, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.25 n.8, p.944-958, August
Gianluca Pignalberi , Rita Cucchiara , Luigi Cinque , Stefano Levialdi, Tuning range image segmentation by genetic algorithm, EURASIP Journal on Applied Signal Processing, v.2003 n.1, p.780-790, January
Geovany Araujo Borges , Marie-Jos Aldon, Line Extraction in 2D Range Images for Mobile Robotics, Journal of Intelligent and Robotic Systems, v.40 n.3, p.267-297, July 2004
Jun Liu , Jim P. Y. Lee , Lingjie Li , Zhi-Quan Luo , K. Max Wong, Online Clustering Algorithms for Radar Emitter Classification, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.27 n.8, p.1185-1196, August 2005
Martin H. C. Law , Mario A. T. Figueiredo , Anil K. Jain, Simultaneous Feature Selection and Clustering Using Mixture Models, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.26 n.9, p.1154-1166, September 2004
Kuo-Lung Wu , Miin-Shen Yang, Alternative learning vector quantization, Pattern Recognition, v.39 n.3, p.351-362, March, 2006
Kuo-Lung Wu , Miin-Shen Yang, Mean shift-based clustering, Pattern Recognition, v.40 n.11, p.3035-3052, November, 2007
Ozy Sjahputera , James M. Keller , J. Wade Davis , Kristen H. Taylor , Farahnaz Rahmatpanah , Huidong Shi , Derek T. Anderson , Samuel N. Blisard , Robert H. Luke , Mihail Popescu , Gerald C. Arthur , Charles W. Caldwell, Relational Analysis of CpG Islands Methylation and Gene Expression in Human Lymphomas Using Possibilistic C-Means Clustering and Modified Cluster Fuzzy Density, IEEE/ACM Transactions on Computational Biology and Bioinformatics (TCBB), v.4 n.2, p.176-189, April 2007
Mohamed Ben Hadj Rhouma , Hichem Frigui, Self-Organization of Pulse-Coupled Oscillators with Application to Clustering, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.23 n.2, p.180-195, February 2001
Jian Yu, General C-Means Clustering Model, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.27 n.8, p.1197-1211, August 2005
Anil K. Jain , Robert P. W. Duin , Jianchang Mao, Statistical Pattern Recognition: A Review, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.22 n.1, p.4-37, January 2000
Jasit S. Suri , Sameer Singh , S. K. Setarehdan , Rakesh Sharma , Keir Bovis , Dorin Comaniciu , Laura Reden, A note on future research in segmentation techniques applied to neurology, cardiology, mammography and pathology, Advanced algorithmic approaches to medical image segmentation: state-of-the-art application in cardiology, neurology, mammography and pathology, Springer-Verlag New York, Inc., New York, NY, 2001
robust clustering;robust statistics;motion estimation;optimal number of clusters;range image segmentation;linear regression;fuzzy clustering;competitive clustering