Mean Shift: A Robust Approach Toward Feature Space Analysis

“Mean Shift: A Robust Approach Toward Feature Space Analysis” by Dorin Comaniciu of Siemens Corporate Research in Princeton and Peter Meer of Rutgers University appeared in IEEE Transactions on Pattern Analysis and Machine Intelligence in May 2002 (vol. 24, no. 5, pp. 603–619), after submission in January 2001. It took a mode-seeking procedure from 1975, set it in a general kernel framework, stated a convergence theorem for kernels with convex profiles, related it to kernel regression and robust M-estimation, and used it in the joint space of pixel position and colour for edge-preserving smoothing and image segmentation, releasing the code. The method as understood today, with its derivation, worked examples and the later convergence results, is in the mean shift article; this article covers the paper as published. Section, equation and figure numbers below come from the authors’ copy on their Rutgers laboratory’s page, an 18-page typeset version with the journal’s running heads, and the corrections from the authors’ note on the same page.

Problem

The paper begins from a complaint about low-level vision: its results depend on the user guessing tuning parameters, while a low-level stage that is to be driven by higher-level information should be reliable and controlled by very few parameters with intuitive meaning (§1). The authors’ answer is feature space analysis: map small subsets of the input (single pixels, or quasi-randomly chosen sets of points as in the probabilistic Hough transform) to points in a parameter space, where significant features become dense regions, that is, clusters. Pooling all the evidence makes the analysis tolerant of noise that would defeat local decisions; the cost is that features with little support in the space can be missed, which the authors propose to offset by adding spatial coordinates to the feature space or by post-processing in the image (§1).

Real feature spaces, the paper argues, defeat most clustering methods. Its Fig. 1 maps a $400 \times 276$ colour photograph to 110,400 points in the $L^*u^*v^*$ colour space, where clusters from the dominant colours blend continuously into each other. Tiling such a space with ellipses introduces severe artifacts; fitting a Gaussian mixture is, in the authors’ words, “doomed to fail” here, even with contaminated (robust) Gaussian components; and mixtures need the number of clusters as an input. Methods that assume a number of clusters or a common cluster shape are dismissed as inadequate. Of the nonparametric alternatives, hierarchical clustering is expensive and lacks a natural stopping rule, which leaves density estimation: dense regions are the modes of the underlying density, and each cluster is delineated around its mode (§1).

For the mode finder, the paper takes the mean shift procedure, which Fukunaga and Hostetler proposed in 1975 and which, in its account, was largely forgotten until Cheng (1995) revived interest. The authors remark that the procedure did not seem to be known in the statistical literature, although Silverman’s book discusses Fukunaga and Hostetler, and that the advantages of a mean shift type procedure in density estimation had only recently been rediscovered (§1).

The paper continues the authors’ own line of work, which it describes in §4 and the acknowledgments. Their first segmentation algorithm (CVPR 1997) clustered colours alone in $L^*u^*v^*$ and had been built into other groups’ systems for image retrieval, face tracking, MPEG-4 object coding, shape detection and texture analysis. Parts of the 2002 material had appeared at ICCV 1999 and CVPR 2000.

Contribution

As the abstract and introduction present them:

• A general nonparametric technique for analysing a multimodal feature space and delineating clusters of arbitrary shape, with the mean shift as its computational module.
• A convergence proof for discrete data: in the abstract’s words, the recursive procedure converges to the nearest stationary point of the underlying density, which makes it usable for detecting modes (see Corrections for what the theorem and its proof support).
• Connections to the Nadaraya–Watson estimator of kernel regression and to robust M-estimators of location.
• Two low-level vision algorithms, discontinuity-preserving smoothing and segmentation, for grey or colour images, in which the only user-set parameter is the resolution of the analysis.

By our reading of the ICCV 1999 text, the conference paper already contained the joint-domain filtering and segmentation and the cameraman, baboon, MIT, lake and hand experiments, with a convergence theorem for the Epanechnikov estimate under uniform-kernel mean shift. The journal version adds the general kernel and profile formulation, Theorem 1 for any convex decreasing profile, the smooth trajectory theorem, the kernel regression and M-estimator sections, the bandwidth-selection techniques of §3.1, the clustering of 110,400 colour points in Fig. 2, and the room, landscape and other segmentation figures.

Method

Kernels, profiles and the mean shift vector

Section 2 builds the procedure from the multivariate kernel density estimate with a bandwidth matrix $\mathbf{H}$ (eqs. 1–2), kernel conditions (3), and two ways of building a multivariate kernel from a univariate one, as a product or by rotation (4). It then restricts itself to radially symmetric kernels $K(\mathbf{x}) = c_{k,d}\, k(\lVert \mathbf{x} \rVert^2)$, defined by their profile $k$ (5), and to a single bandwidth, $\mathbf{H} = h^2 \mathbf{I}$, which requires the feature space to be Euclidean (6). Two profiles are singled out: the Epanechnikov profile, whose kernel minimizes the asymptotic mean integrated squared error (AMISE) of the estimate (7–8), and the normal profile $k_N(x) = \exp(-x/2)$ (9–10), often truncated to finite support.

With $g = -k'$ (13) and the kernel $G$ built from $g$ (14), the gradient of the estimate factors into a density estimate computed with $G$ times the mean shift $\mathbf{m}_{h,G}(\mathbf{x})$, the difference between the $G$-weighted mean of the samples and the window’s centre (15–17). Rearranged, $\mathbf{m}_{h,G}(\mathbf{x}) = \tfrac{1}{2} h^2 c\, \nabla \hat f_{h,K}(\mathbf{x}) / \hat f_{h,G}(\mathbf{x})$ (19): the mean shift is proportional to the density gradient normalized by the density. The paper calls (19) a more general formulation of a property first remarked by Fukunaga and Hostetler; it takes the term shadow for $K$ relative to $G$ from Cheng, noting that Cheng used it in a somewhat different setting (§2.1). Steps are long where the density is low and short near a mode, so the paper calls the procedure an adaptive gradient ascent method. The derivation is worked through in the mean shift article.

Convergence and mode detection

The procedure moves the window to the weighted mean, $\mathbf{y}_{j+1} = \sum_i \mathbf{x}_i\, g(\lVert (\mathbf{y}_j - \mathbf{x}_i)/h \rVert^2) / \sum_i g(\lVert (\mathbf{y}_j - \mathbf{x}_i)/h \rVert^2)$, eq. (20) as corrected by the authors (see Limitations). Theorem 1 (§2.2) states: if the kernel $K$ has a convex and monotonically decreasing profile, the sequences $\{\mathbf{y}_j\}$ and $\{\hat f_{h,K}(\mathbf{y}_j)\}$ converge, and the density values increase monotonically. The paper says the theorem generalizes the ICCV 1999 result, in which $K$ was the Epanechnikov kernel and $G$ the uniform kernel, that it holds when each sample carries a nonnegative weight, and that the first author’s thesis gives an example of nonconvergence for a nonconvex $K$. It distinguishes its procedure from Cheng’s: in its account, almost all of Cheng’s convergence discussion concerns the blurring process, in which the data themselves move at each step, and Cheng attributed the convergence of the non-blurring procedure to the gradient ascent nature of (19). Comaniciu and Meer counter that moving along the gradient guarantees convergence only for infinitesimal steps, and credit the guarantee instead to the adaptive magnitude of the mean shift vector, which removes the need for a step-size rule. With a uniform $G$, convergence takes finitely many steps because only finitely many distinct means exist; with a weighting $G$, the procedure converges only in the limit and is stopped by a lower bound on the step.

The proof (Appendix) uses the convexity of the profile to show that each step raises the density estimate by at least a positive multiple of the squared step length (A.1–A.6), and then sums these bounds over $m$ consecutive steps to conclude that $\{\mathbf{y}_j\}$ is a Cauchy sequence (A.7).

From Theorem 1, §2.3 concludes that the mean shift vector tends to zero, so the limit $\mathbf{y}_c$ is a stationary point of the estimate (22–23), and that the monotone increase of the density satisfies the conditions of the Capture Theorem from Bertsekas’s Nonlinear Programming, so trajectories are attracted by local maxima that are unique stationary points in a small neighbourhood. Mode detection therefore runs the procedure to find stationary points and keeps only local maxima, tested by perturbing each stationary point by a small random vector and checking that the procedure returns to it. The set of points converging to a mode is its basin of attraction.

Theorem 2 (§2.4) adds that with the normal kernel the cosine between consecutive mean shift vectors is strictly positive (25), so the path to the mode is smooth, unlike the zigzag of steepest ascent with line search on surfaces with narrow valleys. The authors conclude that the normal kernel seems the best choice for mean shift, but report that it needs many steps in practice, so most of their experiments use the uniform kernel, with finite convergence, while noting that results almost always improve with the normal kernel.

Kernel regression and M-estimators

Section 2.5 derives (19) a second way. The Nadaraya–Watson estimator (27), which fits a constant to the data in a moving window, has an asymptotic bias proportional to $h^2 f'(x)/f(x)$ when the regression function is the identity (28–29); with $m(x) = x$ the estimator reduces to the mean shift step, so the procedure, as the paper puts it, turns the inherent bias of zero-order kernel regression to its advantage. Section 2.6 shows that the iteration computing a location M-estimate with weight function $w$ (30–32) is identical to (20) when $w \equiv g$, and that the estimate is the mode of the density estimated with $K$ (33–34); the convexity of $k$ required in Theorem 1 matches the usual requirements on the M-estimator’s objective $\rho$. The authors note that only one earlier paper, on edge-preserving smoothers, had discussed the relation between M-estimators and kernel density estimation.

Feature space analysis

Section 3 turns the procedure into a clustering technique. Procedures are started from points that tessellate the occupied feature space (regular tessellations are not required). Two issues come first: the metric of the feature space and the shape of the kernel. A Mahalanobis metric could be absorbed into $\mathbf{H}$, but the authors prefer to make the space Euclidean and use $\mathbf{H} = h^2 \mathbf{I}$; one bandwidth for the whole space requires that local density variations near a significant mode be smaller than the whole support of a significant mode elsewhere. Convergence points closer than the bandwidth are fused, keeping the one with the highest density, which removes the scatter caused by the stopping threshold and by plateaus; valleys along cuts between modes can be tested for significance; and each mode’s basin of attraction delineates a cluster of arbitrary shape, with majority logic near the boundaries. The authors stress that in vision the image itself usually provides an independent check on uncertain local decisions in feature space.

Fig. 2 illustrates this on the first two coordinates ($L^*$, $u^*$) of the Fig. 1 colour space, 110,400 points: 159 uniform-kernel procedures, with trajectories drawn over the Epanechnikov density estimate, produce seven clusters after pruning, and some trajectories are stopped early on local plateaus. For a quantitative check the paper refers to the authors’ 1999 Pattern Analysis and Applications paper, where on a synthetic bimodal normal distribution the technique reached a classification error similar to that of the optimal Bayesian classifier.

Section 3.1 lists four ways to choose the bandwidth: the AMISE-optimal bandwidth (the multivariate formula is of little practical use, while the plug-in rule is reliable in one dimension); the centre of the largest range of bandwidths giving the same number of clusters; maximizing a cluster validity index; and top-down information from the user or a higher-level module. Locally adaptive bandwidths are deferred to the authors’ variable-bandwidth paper (ICCV 2001), summarized as showing that the local bandwidth should maximize the magnitude of the normalized mean shift vector. Section 3.2 speeds up the computation with a binned estimator on a regular grid, multidimensional range searching, and a fast Euclidean distance criterion from vector quantization.

The joint spatial–range domain

Section 4 represents an image as a lattice of $p$-dimensional pixel values ($p = 1$ for grey, 3 for colour, more for multispectral images) and concatenates each pixel’s two lattice coordinates (the spatial domain) with its value (the range domain) into a vector of dimension $d = p + 2$. Since the two parts have different natures, the kernel is a product of two radially symmetric kernels with bandwidths $h_s$ and $h_r$ (35), and the user sets only $\mathbf{h} = (h_s, h_r)$; according to the paper, an Epanechnikov or truncated normal kernel always gives satisfactory performance. Colour is converted to $L^*u^*v^*$ so that Euclidean distance approximates perceived colour difference; the authors found no clear advantage of $L^*u^*v^*$ over $L^*a^*b^*$ and chose the former for a linear mapping property.

Mean shift filtering (§4.1.1) runs a procedure from every pixel in the joint domain and gives the pixel the range part of its convergence point, keeping its position. The paper sets it against adaptive Wiener filtering, diffusion methods, which lack a stopping criterion and eventually collapse the image to a flat surface, and the bilateral filter of Tomasi and Manduchi (1998), which also weights in both domains but uses a fixed window, whereas the mean shift window moves along the joint density gradient and so uses information beyond the initial window. The authors attribute the smoothing quality of both filters to a twofold weighting similar in spirit to generalized (bounded-influence) M-estimators, which down-weight outliers in the input domain as well.

Mean shift segmentation (§4.2.1) runs the filter and stores the full $d$-dimensional convergence point of every pixel, groups all convergence points closer than $h_s$ in space and $h_r$ in range, which concatenates the basins of attraction of nearby modes, labels each pixel by its group, and optionally eliminates regions of fewer than $M$ pixels. The authors note that the grouping can incorporate a priori information, such as physics-based segmentation, and that hierarchical processing of the region adjacency graph can speed it up. They compare the method with an attraction force field (whose vectors they find similar to the spatial part of the joint-domain mean shift vector, but which looks for boundaries where the vectors diverge, while mean shift moves away from them), edge flow, the Gaussian “blob” representation, and a hierarchical nonparametric clustering method that, as they report it, takes several minutes in MATLAB for a 2,000-pixel subsample.

Results

The experiments are qualitative figures; the paper has no tables and no quantitative comparison with other segmenters. All segmentation experiments use uniform kernels.

Figure	Image	Task	Parameters	Reported outcome
3	cameraman, $256 \times 256$ grey	filtering, uniform kernel	$(h_s, h_r) = (8, 4)$	grass almost completely smoothed, tripod and buildings preserved; fractions of a second; 3.06 iterations per pixel on average
4	$40 \times 20$ window of cameraman	filtering and segmentation	$(8, 4)$ for the filtering (Fig. 4c)	paths of every other pixel; convergence points in the centre of the plateau, paths on a thin line stay on it
5	baboon, $512 \times 512$ colour	filtering, normal kernels	$h_s \in \{8, 16, 32\}$, $h_r \in \{4, 8, 16\}$	fur texture removed, eyes and whiskers crisp
6	MIT, $256 \times 256$ grey	segmentation	$(h_s, h_r, M) = (8, 7, 20)$	225 regions in fractions of a second
7	room, $256 \times 256$ colour	segmentation	$(8, 5, 20)$	carpet with an illumination gradient kept as one region
8	lake, $512 \times 512$ colour	segmentation	$(16, 7, 40)$	compared with a multiscale method
9	hand, colour	segmentation	$(16, 19, 40)$	textured background handled by large $h_r$ and $M$
10	four landscape images	segmentation	$(8, 7, 100)$	boundaries drawn over the input
11	four other colour images	segmentation	$(8, 7, 20)$	original and segmented side by side
12	cameraman	segmentation and editing	$(8, 4, 10)$	sky and grass replaced with white; supervised texture insertion

The timings are for an optimized C++ implementation on a 600 MHz Pentium III (§4.1); the ICCV 1999 version had reported a few seconds for the same filtering of cameraman with a Java implementation on a 233 MHz Pentium II. The Fig. 5 grid shows eight of the nine bandwidth pairs next to the original, and the text draws the lesson that the spatial bandwidth has a distinct effect from the range bandwidth: with large $h_s$ only features of large spatial support survive, with large $h_r$ only features of high colour contrast, behaviour the authors note was also reported for the bilateral filter.

For segmentation, Fig. 4d shows the grouping step fusing the filtered plateau into larger regions, and the paper says this step adds no significant overhead to the filtering. The comparisons are visual and point to figures in other papers: one-dimensional grey-level clustering and a Gibbs random field method for MIT, recursive thresholding for room, a multiscale edge-and-region method for lake, and region competition and geodesic active contours for hand. The authors’ reading is that joint-domain segmentation keeps regions with a small illumination gradient, such as the sky in MIT and the carpet in room, as single regions, where grey-level or colour clustering alone typically oversegments them.

On parameters, the paper states that the segmentation is not very sensitive to $h_s$ and $h_r$, and that all $256 \times 256$ images used $h_s = 8$ (a $17 \times 17$ window) and all $512 \times 512$ images $h_s = 16$, described as a $31 \times 31$ window (by the $2h_s + 1$ rule that gives 17 for $h_s = 8$, it would be $33 \times 33$; our arithmetic). The range bandwidth $h_r$ and the minimum region size $M$ control the number of regions, and must grow as the image departs from the piecewise-constant model: the heavily textured background of hand needed $h_r = 19$ and $M = 40$, against $h_r = 5$ and $M = 20$ for room. Fig. 12 suggests image editing as an application. The code for both algorithms, integrated into one system with a graphical interface, was released through the laboratory’s web page; the laboratory’s page names it the Edge Detection and Image SegmentatiON (EDISON) system, which also contains an edge detector from a separate paper.

Impact

Carreira-Perpiñán (2007, §1) credits Comaniciu and Meer with a convergence proof for certain isotropic kernels and with applying mean shift to filtering and segmentation, which, in his account, set off a wave of interest in the algorithm for these and other vision applications such as tracking. The tracking work that the paper describes in §5, published at CVPR 2000 by Comaniciu, Ramesh and Meer and listed by the 2002 paper among its preliminary versions, is listed by the PAMI Technical Committee among the 2010 winners of the Longuet-Higgins Prize, which recognizes CVPR papers from ten years earlier with significant impact on computer vision research; the prize belongs to that CVPR paper, not to the 2002 journal paper.

The segmentation algorithm became a standard baseline. In the benchmark of Arbeláez, Maire, Fowlkes and Malik (2011), “Mean Shift” refers to this paper. On the BSDS500 boundary benchmark it reached an F-measure of 0.64 at the optimal dataset scale (ODS), 0.68 at the optimal image scale (OIS) and an average precision of 0.56, level on ODS with the multiscale normalized cuts baseline and ahead of Felzenszwalb–Huttenlocher (0.61), against 0.73, 0.76 and 0.73 for the authors’ gPb-owt-ucm and 0.80 for humans (their Table 1). On the BSDS500 region benchmarks its segment covering at ODS was 0.54, against 0.59 for gPb-owt-ucm and 0.52 for Felzenszwalb–Huttenlocher (Table 2). In their pairwise per-image comparison on BSDS300, all differences between methods were statistically significant except the one between mean shift and normalized cuts (Fig. 18). The normalized cuts paper article reports the same benchmark from the other side.

Mean shift is also listed among the mode-seeking methods for superpixels, and the procedure, beyond segmentation, runs through object tracking, detection merging and embedding clustering, as described in the mean shift article.

Limitations

The paper states several limitations itself:

• Bandwidth. The resolution parameter is left to the user or to a higher level; one global bandwidth requires local density variations to be smaller than the support of other modes (§3), and the authors leave the combination of automatic bandwidth selection with top-down information as an open problem (§5).
• Plateaus and slow kernels. Trajectories can stop prematurely on flat regions (Fig. 2), and the normal kernel, though it gives better results, needs many iterations (§2.4).
• Dimension. Kernel density estimation does not scale well with dimension because of the empty space phenomenon, so feature spaces with more than about six dimensions should be approached carefully, possibly with projection pursuit (§5).
• Speed of segmentation. Unlike the tracker, which ran at 30 frames per second, it was not clear whether segmentation could be made as fast (§5).
• Evaluation. The authors write that the quality of a segmentation can be assessed only in the context of the whole vision task, and choose parameters by that criterion; the paper reports no quantitative segmentation results.

Corrections

The authors’ correction note, linked as errata from their laboratory’s page, states that equations (20) and (24) were misprinted in the published paper. In the authors’ copy, (20) has $\mathbf{x}$ where $\mathbf{y}_j$ belongs inside $g$, and (24), the mean shift step for the normal kernel, has $\mathbf{x}$ for $\mathbf{y}_j$ and no minus sign in the exponent; the note restores both. The corrected (24) uses $\exp(-\lVert (\mathbf{y}_j - \mathbf{x}_i)/h \rVert^2)$, as does the proof of Theorem 2 (B.1), without the factor $\tfrac{1}{2}$ of the normal profile (9); in (24), therefore, $h$ equals $\sqrt{2}$ times the Gaussian’s standard deviation rather than the standard deviation itself (our reading of the two formulas).

The proof of Theorem 1 has a gap that the authors did not correct. Steps (A.1)–(A.6) establish that the density values increase and converge and that the steps shrink, but (A.7) bounds $\lVert \mathbf{y}_{j+m} - \mathbf{y}_j \rVert^2$ by the sum of the squared steps, which would need a triangle inequality for squared distances. Aliyari Ghassabeh (2013) points out that the sequence of points need not be Cauchy and that the proof’s crucial step is therefore incorrect, and the errata of Arias-Castro, Mason and Pelletier (2016), correcting the same mistake in their own paper, name Comaniciu and Meer (2002) among those who made it. As Aliyari Ghassabeh describes the later work, Li, Hu and Wu (2007) proved convergence assuming finitely many stationary points of the estimate in the convex hull of the data, an assumption not shown for common kernels; Carreira-Perpiñán’s EM interpretation does not settle the question either, because EM iterations need not converge in general; and Aliyari Ghassabeh himself proved convergence in one dimension for convex, differentiable, strictly decreasing profiles whose $g$ is also strictly decreasing (his Theorem 1). The mean shift article gives the corrected account.

The abstract’s phrase “the nearest stationary point” is also stronger than the theorem, which shows convergence to a stationary point. Carreira-Perpiñán (2007) found that the convergence domains of Gaussian mean shift can be nonconvex and disconnected, with fractal behaviour, so the mode reached from a starting point is not necessarily the nearest one (our inference from his result).

What Came After

The authors themselves extended the paper in two directions it names: variable-bandwidth mean shift with data-driven scale selection (ICCV 2001) and kernel-based tracking (CVPR 2000, with a journal version in 2003). On the theory side, Carreira-Perpiñán (2007) showed that Gaussian mean shift is an expectation–maximization algorithm and analysed its convergence rate, typically linear and often very slow, and the convergence of the point sequence itself was still being settled in the 2010s. In image segmentation, benchmarks with human ground truth replaced the side-by-side figures of the paper, and contour-based hierarchical methods such as gPb-owt-ucm overtook mean shift on them. Mean shift survives as a building block wherever a density’s modes are wanted without fixing their number: colour clustering, superpixels, merging of detections, and grouping of learned embeddings.