Mean Shift

Mean shift finds the modes, the local maxima, of a probability density known only through samples. It places a window at a starting point, replaces the point by the weighted mean of the samples inside the window, and repeats until the point stops moving. Each step follows the gradient of a kernel density estimate, with a step length that adapts to the density, so the iteration climbs to the nearest peak without a step size to tune. Run from every sample, it groups the samples by the peak they reach, which makes it a clustering method that needs no number of clusters, only a window size, the bandwidth. In computer vision the samples are usually pixels: pixel colours for colour clustering, pixel position and colour together for edge-preserving smoothing and image segmentation, and pixel positions weighted by how well their colour fits a target for object tracking.

The procedure goes back to Fukunaga and Hostetler (1975), was revived by Cheng (1995), and was made a general tool for vision by Comaniciu and Meer (2002), whose formulation this article follows. It builds on probability densities and on gradient ascent.

Problem

Given $n$ samples $\mathbf{x}_1, \dots, \mathbf{x}_n \in \mathbb{R}^d$, stored as the rows of $X \in \mathbb{R}^{n \times d}$, find the modes of the density they were drawn from, and assign each sample to a mode. The number of modes is not known, and the clusters around them may have any shape. Comaniciu and Meer (2002, §1) argue that real feature spaces, such as the colours of a photograph in a perceptually uniform colour space, are of this kind: clusters arise from the dominant colours and blend continuously into each other, so a decomposition into elliptical pieces fits them poorly, and methods that fix the number of clusters or assume one cluster shape are inadequate.

Inputs and Outputs

• Input: the samples $X$; a kernel (in practice flat or Gaussian); a bandwidth $h$, or one bandwidth per group of features; starting points, usually every sample or a coarse grid of them; a stopping tolerance.
• Output: the modes found, each with its basin of attraction, the set of starting points that converge to it. For clustering, a label per sample. For filtering, a new value per pixel; for segmentation, a region label per pixel; for tracking, the target’s position in each frame.

Intuition

Put a window of fixed radius anywhere among the samples. The mean of the samples inside it lies, on average, toward the side where the samples are denser, because that side contributes more of them. Moving the window to that mean and repeating carries it uphill in density until the samples inside are balanced around its centre: a peak. In sparse regions the few samples inside the window pull its centre far, so steps are long where nothing of interest happens; near a peak the samples surround the centre and steps become short. The vector from the window’s centre to the mean of its samples is the mean shift.

Every sample that climbs to the same peak belongs to the same cluster. Nothing in the procedure fixes how many peaks there are; the bandwidth does it implicitly: a narrow window sees every noise bump as a peak, a wide one smooths neighbouring peaks into one.

Algorithm

Mode seeking from one starting point $\mathbf{y}_1$, with a kernel profile $g$ described below:

1. Compute the weighted mean of the samples around the current point, $\mathbf{y}_{j+1} = \sum_i \mathbf{x}_i\, g(\lVert (\mathbf{y}_j - \mathbf{x}_i)/h \rVert^2) \big/ \sum_i g(\lVert (\mathbf{y}_j - \mathbf{x}_i)/h \rVert^2)$.
2. If $\lVert \mathbf{y}_{j+1} - \mathbf{y}_j \rVert$ is below a tolerance, stop: $\mathbf{y}_{j+1}$ is taken as a stationary point of the density estimate.
3. Otherwise set $j \leftarrow j + 1$ and go to step 1.

Clustering by mean shift (Comaniciu & Meer, 2002, §2.3, §3):

1. Run the procedure from starting points that cover the occupied part of the feature space, for example every sample.
2. Fuse convergence points closer than the bandwidth, keeping the one with the highest density, since iterations stopped by a tolerance or on a flat plateau end at slightly different places.
3. Optionally keep only local maxima: perturb each convergence point slightly and rerun; a maximum attracts the iteration back, a saddle point does not.
4. Label each sample by the mode its own iteration reached.

The data never move. A different procedure, which Cheng (1995) called the blurring process, replaces every sample by its mean shift at each iteration, so the data set itself contracts (Comaniciu & Meer, 2002, §2.2; Carreira-Perpiñán, 2007, §1).

Mathematical Formulation

Kernel density estimation

A kernel density estimate, also called a Parzen window estimate, places a bump of fixed width on every sample and adds them (Bishop, 2006, §2.5.1). With a radially symmetric kernel written through its profile $k$, a function on $[0, \infty)$ with $K(\mathbf{x}) = c_{k,d}\, k(\lVert \mathbf{x} \rVert^2)$, and one bandwidth $h$ for all dimensions, the estimate is (Comaniciu & Meer, 2002, eq. 11)

$$\hat p_K(\mathbf{x}) = \frac{c_{k,d}}{n h^d} \sum_{i=1}^{n} k\!\left( \left\lVert \frac{\mathbf{x} - \mathbf{x}_i}{h} \right\rVert^2 \right),$$

where $c_{k,d}$ makes the kernel integrate to one. Two profiles cover most uses. The Epanechnikov profile $k_E(x) = 1 - x$ for $0 \le x \le 1$ and $0$ beyond gives the kernel that minimizes the asymptotic mean integrated squared error of the estimate; the normal profile $k_N(x) = \exp(-x/2)$ gives the Gaussian kernel (Comaniciu & Meer, 2002, eqs. 7–10). The bandwidth sets the smoothness: Bishop (§2.5.1) shows the same data with a small $h$ giving a very noisy estimate and a large $h$ washing out the true bimodality. The single bandwidth also assumes the features are measured on a common Euclidean scale; a full bandwidth matrix could absorb a Mahalanobis metric, but Comaniciu and Meer recommend making the feature space Euclidean instead (§2, §3).

The mean shift is a normalized gradient

The gradient of the estimate, with $g(x) = -k'(x)$ and $g_i = g(\lVert (\mathbf{x} - \mathbf{x}_i)/h \rVert^2)$, factors into two terms (Comaniciu & Meer, 2002, eqs. 12–17):

$$\nabla \hat p_K(\mathbf{x}) = \frac{2 c_{k,d}}{n h^{d+2}} \left[ \sum_{i=1}^{n} g_i \right] \left[ \frac{\sum_{i=1}^{n} \mathbf{x}_i\, g_i}{\sum_{i=1}^{n} g_i} - \mathbf{x} \right].$$

The first bracket is proportional to a density estimate $\hat p_G$ computed with the kernel $G(\mathbf{x}) = c_{g,d}\, g(\lVert \mathbf{x} \rVert^2)$; the second is the mean shift $\mathbf{m}_{h,G}(\mathbf{x})$, the difference between the $G$-weighted mean of the samples and the window’s centre. Rearranged (their eq. 19),

$$\mathbf{m}_{h,G}(\mathbf{x}) = \frac{1}{2} h^2 c\, \frac{\nabla \hat p_K(\mathbf{x})}{\hat p_G(\mathbf{x})}, \qquad c = \frac{c_{g,d}}{c_{k,d}}.$$

The mean shift computed with $G$ points along the gradient of the density estimate built with $K$, divided by a density estimate. $K$ is called the shadow of $G$, a term Comaniciu and Meer take from Cheng (1995). The Epanechnikov kernel is the shadow of the flat (uniform) kernel, since $-k_E'$ is constant on $[0, 1]$: averaging the samples inside a ball climbs the Epanechnikov density estimate. The normal kernel is its own shadow, because the derivative of $\exp(-x/2)$ is proportional to itself; then $\hat p_G = \hat p_K$, $c = 2$, and

$$\mathbf{m}_{h,N}(\mathbf{x}) = h^2\, \nabla \log \hat p_K(\mathbf{x}),$$

a relation the first code example below checks numerically.

An adaptive gradient ascent

The iteration $\mathbf{y}_{j+1} = \mathbf{y}_j + \mathbf{m}_{h,G}(\mathbf{y}_j)$ is therefore gradient ascent on $\hat p_K$ with step length $\tfrac{1}{2} h^2 c / \hat p_G(\mathbf{y}_j)$, large where the density is low and small near a peak; Comaniciu and Meer (§2.1) call it an adaptive gradient ascent method. For the Gaussian kernel it reads $\mathbf{y}_{j+1} = \mathbf{y}_j + h^2 \nabla \log \hat p_K(\mathbf{y}_j)$, gradient ascent on the log-density with the fixed step $h^2$ (this follows from the relation above). Fixed-step gradient methods can overshoot and diverge when the step is too large, which is why gradient descent needs a step-size rule or a line search. Mean shift needs neither: Comaniciu and Meer (§2.2) attribute its guaranteed convergence to the adaptive magnitude of the step, not to the gradient direction alone.

Convergence

Comaniciu and Meer (2002, §2.2, Theorem 1) state the following: if the kernel $K$ has a convex and monotonically decreasing profile, the sequence of points $\{\mathbf{y}_j\}$ and the sequence of density values $\{\hat p_K(\mathbf{y}_j)\}$ converge, and the density values increase monotonically; the Epanechnikov and normal profiles both qualify. Their proof (Appendix) uses the convexity of the profile to show that each step raises the density estimate by at least a positive multiple of $\lVert \mathbf{y}_{j+1} - \mathbf{y}_j \rVert^2$, so the bounded, increasing density values converge and the steps shrink to zero. The last step of the proof, from there to the convergence of the points themselves, does not hold in general: it 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 a Cauchy sequence, and the errata of Arias-Castro, Mason and Pelletier (2016) name this mistake in Comaniciu and Meer’s proof. Convergence of the points has been proved under additional conditions, for example when the estimate has finitely many stationary points (Li, Hu and Wu, 2007, as Aliyari Ghassabeh describes them) or in one dimension (Aliyari Ghassabeh, 2013); in practice the iteration is stopped by a tolerance on the step in any case. With the flat kernel $G$ the procedure stops after finitely many steps, because only finitely many subsets of the samples can fall inside the window; with a kernel that weights samples by distance, such as the Gaussian, it converges only in the limit, and is stopped when the step falls below a tolerance (§2.2). With the normal kernel, consecutive mean shift vectors always make an angle of less than 90° (their Theorem 2), so the path to the mode is smooth, unlike the zigzag of steepest ascent with line search.

Where the points converge, the limit is a stationary point of $\hat p_K$, where $\nabla \hat p_K = \mathbf{0}$; that is why step 3 of the clustering algorithm checks for maxima. Carreira-Perpiñán (2007) showed that mean shift with a Gaussian kernel is an expectation–maximization algorithm, and with other kernels a generalized EM algorithm, also for kernels with a different full covariance per sample, and concluded that it converges from almost any starting point (§1–§3); Aliyari Ghassabeh (2013) notes that EM iterations need not converge in general, so this does not settle the convergence of the points either. In principle it can converge to a saddle point or a minimum, but these are unstable and a small perturbation moves the iteration off them (§2). The EM view also explains why mean shift can be slow: convergence is in general linear, and for Gaussian mean shift the rate is often close to 1 for intermediate bandwidths and the convergence sublinear at bandwidths where two modes merge, a situation that occurs in image segmentation (§4).

Mean Shift Clustering

As a clustering method, mean shift takes the modes of the density as cluster centres and the basins of attraction as clusters, which can have any shape (Comaniciu & Meer, 2002, §3). It needs no number of clusters, has no random initialization, and returns the same result on every run. The bandwidth replaces the number of clusters: Comaniciu and Meer (§3.1) list four ways to choose it: a statistical rule that balances the bias and variance of the density estimate (the plug-in rule, reliable in one dimension, while the multivariate formula is of little practical use); the centre of the largest range of bandwidths over which the number of clusters stays the same; the maximum of a cluster validity index; and the needs of the task. Using one bandwidth for the whole space requires that local density variations near one significant mode be smaller than the whole support of another (§3); they treat locally adaptive bandwidths in a separate paper.

Mean shift sits among the density-based clustering methods. k-means fixes the number of clusters and assigns each sample to the nearest of $k$ centres, so its clusters are the cells of a Voronoi partition; a Gaussian mixture model fixes the number of components and their elliptical shape. Mean shift fixes neither, at the price of the bandwidth and of more computation.

scikit-learn’s MeanShift

In scikit-learn 1.7.2, MeanShift uses a flat kernel, so each step moves a point to the plain mean of the samples within distance bandwidth; the defaults are bandwidth=None, seeds=None, bin_seeding=False, min_bin_freq=1, cluster_all=True and max_iter=300 (scikit-learn documentation). In the source of that version:

• Seeds. By default every sample is a seed. With bin_seeding=True, the seeds are the centres of a grid of cell size bandwidth, keeping cells with at least min_bin_freq samples, which is faster.
• Stopping and merging. A seed stops when its step is at most $10^{-3}$ times the bandwidth or after max_iter iterations. Converged points are sorted by the number of samples within the bandwidth, and any point closer than the bandwidth to a better-supported one is removed.
• Labels. Each sample gets the label of the nearest centre, not of the mode its own iteration reached. With cluster_all=False, samples farther than the bandwidth from every centre are labelled −1. The user guide states that new samples are labelled by the nearest centre as well.
• Bandwidth. With bandwidth=None, estimate_bandwidth(X) is called with quantile=0.3. The function computes, for every sample, the distance to its $\lfloor 0.3\, n \rfloor$-th nearest neighbour (counting the sample itself) and returns the mean of these distances, at least quadratic time in $n$. Its docstring describes quantile=0.5 as using the median of all pairwise distances, which is the per-sample median distance averaged over samples rather than the median over all pairs. The function’s signature sets random_state=0, while its docstring says the default is None; the random state is used only when n_samples subsamples the data.

Filtering and Segmentation in the Joint Domain

An image becomes a set of samples in the joint spatial–range domain by concatenating each pixel’s position $\mathbf{x}^s$ (2 coordinates) with its value $\mathbf{x}^r$ (1 for grey level, 3 for colour), so $d = p + 2$. Position and colour have different units, so Comaniciu and Meer (2002, §4, eq. 35) use the product of two radially symmetric kernels with their own bandwidths,

$$K_{h_s, h_r}(\mathbf{x}) = \frac{C}{h_s^2 h_r^p}\, k\!\left( \left\lVert \frac{\mathbf{x}^s}{h_s} \right\rVert^2 \right) k\!\left( \left\lVert \frac{\mathbf{x}^r}{h_r} \right\rVert^2 \right),$$

where $h_s$ is the spatial bandwidth in pixels and $h_r$ the range bandwidth in colour units. They work in the $L^*u^*v^*$ space for colour images, so that Euclidean distances approximate perceived colour differences (§4; see colour spaces). With a flat kernel, a pixel is in the window when it lies within $h_s$ in position and within $h_r$ in colour.

Mean shift filtering (§4.1.1) runs the procedure from every pixel in the joint domain and gives the pixel the range part of the point it converges to, keeping its position. A window that starts on one side of an edge contains, in the joint domain, only the pixels on that side, because pixels across the edge differ by more than $h_r$; it therefore climbs to the mode of its own region and never averages across the edge. This makes it a discontinuity-preserving smoother, related to the bilateral filter, which weights pixels by spatial and range distance in the same way but averages over a window fixed at the pixel; the mean shift window moves, so information from beyond the initial window is used (§4.1).

Mean shift segmentation (§4.2.1) groups the pixels whose convergence points lie closer than $h_s$ in space and $h_r$ in range, which joins the basins of attraction of nearby modes, labels each pixel by its group, and optionally removes regions of fewer than $M$ pixels. The three parameters $(h_s, h_r, M)$ are the whole user interface: $h_r$ and $M$ control the number of regions, and the more an image departs from piecewise-constant colour, as in textured regions, the larger they must be (Comaniciu & Meer, 2002, §4.2). Among methods for superpixels, mean shift is a mode-seeking method that does not control the number or compactness of its regions.

Mean Shift Tracking

Comaniciu, Ramesh and Meer (2000; 2003) turned mean shift into a tracker by making the density to be climbed a measure of how well a window’s contents match a target. The target model is an $m$-bin colour histogram $\hat{\mathbf{q}}$ of the pixels in an elliptical region, normalized to a unit circle, with each pixel weighted by a kernel profile $k$ that gives less weight to pixels near the boundary, where occlusion and background are most likely (2003, §2.1). A candidate at position $\mathbf{y}$ in a new frame has the histogram $\hat{\mathbf{p}}(\mathbf{y})$, computed the same way from the pixels $\mathbf{x}_i$ within bandwidth $h$ of $\mathbf{y}$ (§2.2). Similarity is the Bhattacharyya coefficient (§3),

$$\rho(\mathbf{y}) = \sum_{u=1}^{m} \sqrt{\hat p_u(\mathbf{y})\, \hat q_u},$$

the cosine of the angle between the vectors of square-rooted histograms, with $\sqrt{1 - \rho}$ as the distance (see similarity and distance measures). Because the kernel weighting makes $\rho$ a smooth function of $\mathbf{y}$, it can be maximized by gradient ascent. Expanding $\rho$ to first order around the previous position $\hat{\mathbf{y}}_0$ leaves, up to a constant, a kernel density estimate over pixel positions in which pixel $i$, falling in histogram bin $b(\mathbf{x}_i)$, carries the weight (§4.1, eqs. 9–10)

$$w_i = \sqrt{\frac{\hat q_{b(\mathbf{x}_i)}}{\hat p_{b(\mathbf{x}_i)}(\hat{\mathbf{y}}_0)}} .$$

Pixels whose colour is more common in the target than in the current window pull the window toward them. Mean shift on this weighted density moves the window to $\hat{\mathbf{y}}_1 = \sum_i \mathbf{x}_i w_i g_i / \sum_i w_i g_i$; with the Epanechnikov profile $g$ is constant, and the new position is simply the $w$-weighted mean of the pixel positions (eqs. 11–13). The tracker iterates from the previous frame’s position until the shift falls within a pixel, with a cap on the number of iterations, typically 20 (§4.2). The algorithm includes a safeguard that halves the move when $\rho$ decreases, but the authors found the coefficient failed to increase in only 0.1% of cases and drop the step in practice (§4.1–§4.2). Scale is adapted by rerunning the search with the bandwidth changed by a typical ±10% and smoothing the best bandwidth over time (§4.3). The window can find the target as long as the target’s centre stays within the area the target covered in the previous frame, the “operational basin of attraction” (§4.2).

CamShift, the Continuously Adaptive Mean Shift of Bradski (1998), is an earlier mean shift tracker for faces. Each frame is turned into a probability image by looking up every pixel’s hue in a skin-colour hue histogram; flat-kernel mean shift then moves a search window to the centroid of the probability mass inside it, and the window’s size is reset from the zeroth moment, the total probability under it, so the window grows and shrinks with the face. The second moments give the face’s orientation (head roll). Comaniciu, Ramesh and Meer (2003, §7) note that projecting a model histogram directly onto the frame, as Bradski’s face tracker does, can introduce a large bias in the estimated location and gives a scale-variant measure.

Parameters

• Bandwidth $h$ (or $h_s$ and $h_r$). The only essential parameter: it is the scale at which a bump counts as a mode. Too small and noise creates many modes; too large and distinct modes merge. A useful practical check, the second of Comaniciu and Meer’s four rules, is to pick a value inside a range of bandwidths over which the number of modes does not change.
• Kernel profile. The flat kernel (the Epanechnikov shadow) stops in finitely many steps and is cheap; the Gaussian kernel converges only in the limit but follows a smooth path and gives a smooth density estimate.
• Seeds. Every sample, a grid, or a subsample. Fewer seeds save time but can miss small modes, and a seed placed between clusters can end at a point that is not a mode of any of them (shown in the first example).
• Merge radius and tolerance. Converged points closer than the bandwidth are normally fused; the stopping tolerance is a small fraction of the bandwidth.
• Segmentation: the minimum region size $M$. Tracking: the histogram’s bins, the window’s shape and size, and the maximum number of iterations per frame.

Complexity

Each iteration from one starting point visits the samples inside the window; without a spatial index, it costs $O(nd)$ to compute the distances to all $n$ samples, so running from every sample costs $O(n^2 d)$ per iteration (our reasoning). scikit-learn’s flat-kernel implementation finds the neighbours with a ball tree and its documentation gives a cost tending towards $O(T n \log n)$ in low dimensions and $O(T n^2)$ in high dimensions, with $n$ the number of samples and $T$ described as “the number of points”, presumably the seeds (our reading). Comaniciu and Meer (2002, §3.2) sped up the procedure by binning the data on a regular grid and by range searching. In images the joint domain is cheap to search because the spatial neighbours of a pixel are known from the lattice: only the $(2h_s + 1)^2$ pixels around the window need testing.

Kernel density estimation scales poorly with the dimension, mostly because of the empty space phenomenon, by which most of the mass in a high-dimensional space lies in a small region; Comaniciu and Meer (§5) advise care beyond about six dimensions. The tracker is cheap: the cost per frame is a few weighted histograms of the target window, about 24 times less than an exhaustive search of the same area in their analysis for a 50×50-pixel target and 4.19 iterations per frame (Comaniciu, Ramesh & Meer, 2003, §4.4).

Implementation

The three examples below use synthetic images with known regions, colours and motion, and run with NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2 (OMP_NUM_THREADS=1).

Modes of pixel colours

A 64×64 image has five flat colour regions, two of them dark and light red, plus Gaussian noise with a standard deviation of 0.04 per channel. Each pixel’s RGB colour is a sample in $\mathbb{R}^3$. The script checks the gradient relation for the Gaussian kernel, finds the modes with a NumPy implementation and with scikit-learn, and sweeps the bandwidth.

import numpy as np
from sklearn.cluster import MeanShift, estimate_bandwidth

# A 64x64 RGB image with five flat colour regions, plus Gaussian noise.
TRUE = np.array([[0.85, 0.15, 0.15],   # red background
                 [0.55, 0.10, 0.10],   # darker red square, close to the background colour
                 [0.15, 0.55, 0.20],   # green disk
                 [0.20, 0.25, 0.80],   # blue bar
                 [0.90, 0.85, 0.20]])  # yellow triangle
H = W = 64
r, c = np.mgrid[0:H, 0:W]
label = np.zeros((H, W), int)
label[8:28, 36:58] = 1
label[(r - 42) ** 2 + (c - 18) ** 2 < 12 ** 2] = 2
label[50:58, 30:62] = 3
label[(r > 6) & (r < 30) & (c > 4) & (c < 30) & (c - 4 > (r - 6))] = 4
rng = np.random.default_rng(0)
img = TRUE[label] + rng.normal(0, 0.04, (H, W, 3))
X = img.reshape(-1, 3)                       # one pixel colour per row
print("pixels per region:", np.bincount(label.ravel()))

def mean_shift(X, seeds, h, kernel="flat", tol=1e-3, max_iter=300):
    """Move every seed to a mode of the KDE of X (the data never move)."""
    Y = seeds.astype(float).copy()
    for it in range(max_iter):
        d2 = ((Y[:, None, :] - X[None, :, :]) ** 2).sum(-1)
        if kernel == "flat":                 # g = indicator of the ball: Epanechnikov shadow
            w = (d2 <= h * h).astype(float)
        else:                                # g of the normal profile is again a Gaussian
            w = np.exp(-0.5 * d2 / (h * h))
        Ynew = w @ X / w.sum(1, keepdims=True)
        step = np.linalg.norm(Ynew - Y, axis=1)
        Y = Ynew
        if step.max() <= tol * h:
            break
    return Y, it + 1

def merge_modes(Y, X, h):
    """Fuse converged points closer than h, keeping the one with most points in its window."""
    cand = np.unique(np.round(Y, 6), axis=0)
    support = (((cand[:, None] - X[None]) ** 2).sum(-1) <= h * h).sum(1)
    order = np.argsort(-support, kind="stable")
    keep = []
    for i in order:
        if all(np.linalg.norm(cand[i] - cand[j]) >= h for j in keep):
            keep.append(i)
    return cand[keep]

def match(modes, truth):
    d = np.linalg.norm(modes[:, None] - truth[None], axis=2)
    return d.min(0)                          # distance from each true colour to the nearest mode

# 1. The mean shift vector is the normalized density gradient (Gaussian kernel: m = h^2 grad log f).
h = 0.1
def log_kde(x):
    return np.log(np.exp(-0.5 * ((x - X) ** 2).sum(1) / h ** 2).sum())
for x in [np.array([0.5, 0.3, 0.3]), np.array([0.3, 0.4, 0.5])]:
    w = np.exp(-0.5 * ((x - X) ** 2).sum(1) / h ** 2)
    m = w @ X / w.sum() - x
    eps = 1e-5
    g = np.array([(log_kde(x + eps * e) - log_kde(x - eps * e)) / (2 * eps) for e in np.eye(3)])
    print("mean shift", np.round(m, 5), " h^2 * grad log f", np.round(h * h * g, 5))

# 2. Modes of the colour density, flat and Gaussian kernels, against the true colours.
for kernel, h in [("flat", 0.12), ("gauss", 0.06)]:
    Y, n_it = mean_shift(X, X, h, kernel)
    modes = merge_modes(Y, X, h)
    print(f"{kernel:5s} h={h}: {len(modes)} modes after {n_it} iterations; "
          f"max distance true colour -> nearest mode {match(modes, TRUE).max():.4f}")

# 3. Cross-check against scikit-learn's MeanShift (flat kernel, every pixel a seed).
for h in [0.12, 0.4]:
    ms = MeanShift(bandwidth=h).fit(X)
    Y, _ = mean_shift(X, X, h, "flat")
    ours = merge_modes(Y, X, h)
    d = np.linalg.norm(ms.cluster_centers_[:, None] - ours[None], axis=2)
    # scikit-learn labels a pixel by its nearest centre, not by the mode its own iteration reached
    basin = np.argmin(np.linalg.norm(Y[:, None] - ms.cluster_centers_[None], axis=2), axis=1)
    print(f"h={h}: sklearn {len(ms.cluster_centers_)} centres, ours {len(ours)}; "
          f"largest gap {d.min(1).max():.1e}; pixels labelled by nearest centre "
          f"but climbing to another mode: {int((basin != ms.labels_).sum())}")

# 4. Bandwidth decides how many modes there are.
def errors(labels):
    """Pixels not in their cluster's majority true region."""
    return sum(len(m) - np.bincount(m).max() for k in np.unique(labels)
               for m in [label.ravel()[labels == k]])

h_est = float(estimate_bandwidth(X))
print(f"estimate_bandwidth with its defaults: {h_est:.3f}")
for h in [0.03, 0.06, 0.1, 0.2, h_est, 0.4, 0.6]:
    ms = MeanShift(bandwidth=h).fit(X)
    print(f"h={h:4.2f}: {len(ms.cluster_centers_):3d} modes, "
          f"largest clusters {np.sort(np.bincount(ms.labels_))[::-1][:5].tolist()}, "
          f"pixels outside their cluster's majority region {errors(ms.labels_)}")
ms = MeanShift(bandwidth=0.4, bin_seeding=True).fit(X)
print("h=0.40, bin_seeding=True: distance from each centre to the nearest true colour",
      np.round(np.linalg.norm(ms.cluster_centers_[:, None] - TRUE[None], axis=2).min(1), 3))

pixels per region: [2664  440  437  256  299]
mean shift [ 0.08439 -0.16577 -0.16778]  h^2 * grad log f [ 0.08439 -0.16577 -0.16778]
mean shift [-0.10689  0.00479 -0.0147 ]  h^2 * grad log f [-0.10689  0.00479 -0.0147 ]
flat  h=0.12: 5 modes after 7 iterations; max distance true colour -> nearest mode 0.0041
gauss h=0.06: 5 modes after 10 iterations; max distance true colour -> nearest mode 0.0043
h=0.12: sklearn 5 centres, ours 5; largest gap 5.6e-04; pixels labelled by nearest centre but climbing to another mode: 0
h=0.4: sklearn 4 centres, ours 4; largest gap 6.8e-07; pixels labelled by nearest centre but climbing to another mode: 0
estimate_bandwidth with its defaults: 0.270
h=0.03: 210 modes, largest clusters [657, 607, 555, 153, 103], pixels outside their cluster's majority region 0
h=0.06:   8 modes, largest clusters [2664, 440, 425, 299, 238], pixels outside their cluster's majority region 0
h=0.10:   5 modes, largest clusters [2664, 440, 437, 299, 256], pixels outside their cluster's majority region 0
h=0.20:   5 modes, largest clusters [2663, 441, 437, 299, 256], pixels outside their cluster's majority region 1
h=0.27:   4 modes, largest clusters [3104, 437, 299, 256], pixels outside their cluster's majority region 440
h=0.40:   4 modes, largest clusters [3104, 437, 299, 256], pixels outside their cluster's majority region 440
h=0.60:   1 modes, largest clusters [4096], pixels outside their cluster's majority region 1432
h=0.40, bin_seeding=True: distance from each centre to the nearest true colour [0.044 0.298 0.003 0.002]

The first two lines confirm the normalized-gradient relation: with a Gaussian kernel the mean shift vector equals $h^2 \nabla \log \hat p$ to five decimals, even at points far from any mode. Both kernels recover all five region colours to within 0.005, about a tenth of the noise level, from every pixel as a seed. The flat kernel with $h = 0.12$ and the Gaussian kernel with $h = 0.06$ both find the five modes; a Gaussian with standard deviation $h$ weighs samples out to about $2h$, so the two windows are of similar reach (our reading). The NumPy flat-kernel version and MeanShift agree to within 6×10⁻⁴, a few times the stopping tolerance of $10^{-3} h$, and on this image labelling by nearest centre and by basin of attraction coincide.

The bandwidth sweep shows the scale parameter at work. At $h = 0.03$, smaller than the noise, every noise bump is a mode and the image falls into 210 clusters; at 0.06 a few small spurious clusters remain. From 0.1 to 0.2 the count is stable at five and the clusters match the five regions to within one pixel, the stability range that Comaniciu and Meer suggest using. The default estimate_bandwidth gives 0.27, which merges the two reds (their colours are 0.31 apart) and mislabels the 440 dark-red pixels; at 0.6 everything is one cluster. The last line shows a seeding pitfall: with bin_seeding=True at $h = 0.4$, one grid seed between the dark-red and green clusters converges to a point about 0.3 from both of their colours and farther from the others, a stationary point of the flat-kernel estimate where the window holds pixels of both clusters, and the merge step, which keeps the better-supported of two centres closer than the bandwidth, keeps it and drops the green mode, which is missing from the list.

Joint-domain filtering and segmentation

A 64×64 grey image holds a background at 0.2, a square at 0.5 and a disk at 0.8, with Gaussian noise of standard deviation 0.08. The script runs flat-kernel mean shift filtering in the joint domain, compares it with Gaussian blur away from and near the edges, and segments the result for several range bandwidths. The removal of regions under $M$ pixels merges each small region into the adjacent region with the closest mean, one region at a time, our choice for this example.

import numpy as np
from scipy import ndimage

# A 64x64 grey image: background 0.2, a square at 0.5, a disk at 0.8, noise sigma 0.08.
H = W = 64
r, c = np.mgrid[0:H, 0:W]
truth = np.zeros((H, W), int)
truth[10:34, 8:32] = 1
truth[(r - 40) ** 2 + (c - 42) ** 2 < 14 ** 2] = 2
clean = np.array([0.2, 0.5, 0.8])[truth]
rng = np.random.default_rng(1)
noisy = clean + rng.normal(0, 0.08, (H, W))

def ms_filter(img, hs, hr, tol=0.01, max_iter=100):
    """Mean shift filtering in the joint spatial-range domain with a uniform kernel:
    a pixel is in the window when it is within hs in position AND within hr in grey level."""
    pos = np.stack([r.ravel(), c.ravel()], 1).astype(float)
    val = img.ravel().astype(float)
    ys, yr = pos.copy(), val.copy()                  # one mean shift procedure per pixel
    R = int(np.ceil(hs)) + 1
    off = np.array([(i, j) for i in range(-R, R + 1) for j in range(-R, R + 1)])
    active = np.ones(len(val), bool)
    for it in range(max_iter):
        a = np.flatnonzero(active)
        cand = np.round(ys[a])[:, None, :].astype(int) + off[None]      # nearby lattice sites
        ok = (cand[..., 0] >= 0) & (cand[..., 0] < H) & (cand[..., 1] >= 0) & (cand[..., 1] < W)
        cand = np.clip(cand, 0, [H - 1, W - 1])
        idx = cand[..., 0] * W + cand[..., 1]
        inside = ok & (((pos[idx] - ys[a, None]) ** 2).sum(-1) <= hs * hs) \
                    & (np.abs(val[idx] - yr[a, None]) <= hr)
        n = inside.sum(1)
        new_s = (inside[..., None] * pos[idx]).sum(1) / n[:, None]
        new_r = (inside * val[idx]).sum(1) / n
        step = np.sqrt(((new_s - ys[a]) ** 2).sum(1) / hs ** 2 + (new_r - yr[a]) ** 2 / hr ** 2)
        ys[a], yr[a] = new_s, new_r
        active[a[step < tol]] = False
        if not active.any():
            break
    return yr.reshape(H, W), ys, yr, it + 1

def segment(ys, yr, hs, hr, M):
    """Group convergence points closer than hs in space and hr in range; then merge
    regions of fewer than M pixels into the adjacent region with the closest mean."""
    n = len(yr)
    parent = np.arange(n)
    def find(i):
        while parent[i] != i:
            parent[i] = parent[parent[i]]
            i = parent[i]
        return i
    near = (((ys[:, None] - ys[None]) ** 2).sum(-1) < hs * hs) & (np.abs(yr[:, None] - yr[None]) < hr)
    for i, j in zip(*np.nonzero(np.triu(near, 1))):
        parent[find(i)] = find(j)
    lab = np.array([find(i) for i in range(n)]).reshape(H, W)
    # connected components in the image, so a region is one spatial piece
    regions = np.zeros((H, W), int)
    k = 0
    for v in np.unique(lab):
        comp, m = ndimage.label(lab == v)
        regions[comp > 0] = comp[comp > 0] + k
        k += m
    n_before = len(np.unique(regions))
    while True:
        ids, sizes = np.unique(regions, return_counts=True)
        if sizes.min() >= M:
            break
        small = ids[np.argmin(sizes)]
        mask = regions == small
        ring = ndimage.binary_dilation(mask) & ~mask
        nbrs = np.unique(regions[ring])
        means = np.array([yr.reshape(H, W)[regions == q].mean() for q in nbrs])
        regions[mask] = nbrs[np.argmin(np.abs(means - yr.reshape(H, W)[mask].mean()))]
    return regions, n_before

boundary = ndimage.morphological_gradient(truth, size=2) > 0     # pixels touching another region
dist = ndimage.distance_transform_edt(~boundary)
near, far = dist <= 2, dist >= 5

def report(name, out):
    rmse = lambda m: np.sqrt(np.mean((out - clean)[m] ** 2))
    print(f"{name:28s} RMSE: all {rmse(dist >= 0):.4f}, >= 5 px from an edge {rmse(far):.4f}, "
          f"<= 2 px from an edge {rmse(near):.4f}")

hs, hr, M = 5, 0.15, 20
report("noisy input", noisy)
filtered, ys, yr, n_it = ms_filter(noisy, hs, hr)
report(f"mean shift (hs={hs}, hr={hr})", filtered)
for s in [1.0, 2.0]:
    report(f"Gaussian blur sigma={s}", ndimage.gaussian_filter(noisy, s))
print("mean shift iterations until every pixel stopped:", n_it)

for hr_ in [0.05, 0.1, 0.15, 0.2, 0.3]:
    f, ys_, yr_, _ = ms_filter(noisy, hs, hr_)
    regions, n0 = segment(ys_, yr_, hs, hr_, M)
    ids = np.unique(regions)
    wrong = sum((regions == q).sum() - np.bincount(truth[regions == q]).max() for q in ids)
    print(f"hr={hr_:4.2f}: {n0:4d} regions after grouping, {len(ids)} after merging those under "
          f"M={M} pixels; pixels outside their region's majority true label: {wrong}")

noisy input                  RMSE: all 0.0802, >= 5 px from an edge 0.0795, <= 2 px from an edge 0.0836
mean shift (hs=5, hr=0.15)   RMSE: all 0.0184, >= 5 px from an edge 0.0144, <= 2 px from an edge 0.0341
Gaussian blur sigma=1.0      RMSE: all 0.0495, >= 5 px from an edge 0.0349, <= 2 px from an edge 0.0966
Gaussian blur sigma=2.0      RMSE: all 0.0664, >= 5 px from an edge 0.0412, <= 2 px from an edge 0.1377
mean shift iterations until every pixel stopped: 41
hr=0.05:  215 regions after grouping, 2 after merging those under M=20 pixels; pixels outside their region's majority true label: 573
hr=0.10:   20 regions after grouping, 4 after merging those under M=20 pixels; pixels outside their region's majority true label: 6
hr=0.15:    7 regions after grouping, 3 after merging those under M=20 pixels; pixels outside their region's majority true label: 6
hr=0.20:    2 regions after grouping, 2 after merging those under M=20 pixels; pixels outside their region's majority true label: 573
hr=0.30:    2 regions after grouping, 2 after merging those under M=20 pixels; pixels outside their region's majority true label: 573

Mean shift filtering lowers the error both inside the regions and at the edges, while Gaussian blur trades one for the other: a wider blur removes more noise in flat areas but averages across the edges, and within 2 pixels of an edge both blurs end worse than the noisy input. The mean shift windows stop at the edges because pixels across them are farther than $h_r$ in grey level; the error left near the edges presumably comes from noisy pixels whose own value falls within $h_r$ of the other region (our explanation).

The segmentation sweep shows what $h_r$ does. At $h_r = 0.15$ the grouping leaves 7 regions, the merge of regions under 20 pixels leaves 3, and all but 6 pixels agree with the true regions. With $h_r = 0.10$ more fragments survive the grouping, and the background stays split into two large regions, giving 4 regions. With $h_r = 0.05$, smaller than the noise, filtering barely smooths, the grouping leaves 215 fragments, and merging them one by one under the size rule absorbs the square into the background. At $h_r = 0.2$ and above, the square merges with the background although their levels differ by 0.3. The grouping is transitive, so a few edge pixels whose convergence values lie between the two levels are enough to link the regions (our explanation). The range bandwidth has to sit between the noise and the smallest contrast that should survive.

Tracking a coloured disk

A reddish disk of radius 12 moves over a background of blue, green and grey 8×8 blocks along a known path; a reddish rectangle elsewhere in the frame is a distractor. Each frame adds Gaussian noise. The tracker follows Comaniciu, Ramesh and Meer: an 8×8×8 RGB histogram weighted by the Epanechnikov profile, the weights $w_i$, and the weighted mean as the new position. The speed parameter scales how far the disk moves between frames.

import numpy as np

# Synthetic video: a reddish disk (radius 12) moving over a cluttered background.
H, W, R, T = 120, 160, 12, 40
rng = np.random.default_rng(2)
yy, xx = np.mgrid[0:H, 0:W]
blocks = rng.integers(0, 3, (H // 8 + 1, W // 8 + 1))          # blue, green, grey patches
palette = np.array([[0.2, 0.3, 0.7], [0.25, 0.6, 0.3], [0.5, 0.5, 0.5]])
background = palette[blocks[yy // 8, xx // 8]]
background[30:50, 100:125] = [0.85, 0.35, 0.2]                  # a distractor of similar colour

def trajectory(t, speed):
    w = 0.05 * speed                                             # known path, (row, col)
    return np.array([60 + 35 * np.cos(0.7 * w * t), 80 - 55 * np.cos(w * t)])

def frame(centre, rng):
    img = background.copy()
    img[(yy - centre[0]) ** 2 + (xx - centre[1]) ** 2 <= R ** 2] = [0.9, 0.3, 0.25]
    return np.clip(img + rng.normal(0, 0.05, img.shape), 0, 1)

BINS = 8
def bins(img):                                                   # RGB quantized to 8x8x8
    q = np.minimum((img * BINS).astype(int), BINS - 1)
    return (q[..., 0] * BINS + q[..., 1]) * BINS + q[..., 2]

def window(c, h):
    """Pixels inside the circle of radius h around c, with Epanechnikov weights 1 - r^2."""
    r0, c0 = int(round(c[0])), int(round(c[1]))
    ys, xs = np.mgrid[r0 - h:r0 + h + 1, c0 - h:c0 + h + 1]
    d2 = ((ys - c[0]) ** 2 + (xs - c[1]) ** 2) / h ** 2
    keep = (d2 < 1) & (ys >= 0) & (ys < H) & (xs >= 0) & (xs < W)
    return ys[keep], xs[keep], 1 - d2[keep]

def histogram(b, c, h):
    ys, xs, k = window(c, h)
    p = np.bincount(b[ys, xs], weights=k, minlength=BINS ** 3)
    return p / p.sum()

def track(speed, h=R, eps=0.5, max_iter=20):
    rng = np.random.default_rng(3)
    c = trajectory(0, speed)
    q = histogram(bins(frame(c, rng)), c, h)                    # target model from frame 0
    errors, iters = [], []
    for t in range(1, T):
        b = bins(frame(trajectory(t, speed), rng))
        for it in range(1, max_iter + 1):
            p = histogram(b, c, h)
            ys, xs, _ = window(c, h)
            u = b[ys, xs]
            w = np.sqrt(q[u] / np.maximum(p[u], 1e-12))          # weights w_i = sqrt(q_u / p_u(y0))
            if w.sum() == 0:                                     # no target colour in the window
                break
            new = np.array([ys @ w, xs @ w]) / w.sum()           # Epanechnikov profile: g constant
            done = np.linalg.norm(new - c) < eps
            c = new
            if done:
                break
        rho = np.sqrt(histogram(b, c, h) * q).sum()             # Bhattacharyya coefficient
        errors.append(np.linalg.norm(c - trajectory(t, speed)))
        iters.append(it)
    return np.array(errors), np.array(iters), rho, c

for speed in [1.0, 3.0, 5.0, 7.0, 9.0]:
    step = np.linalg.norm(np.diff([trajectory(t, speed) for t in range(T)], axis=0), axis=1).max()
    e, n, rho, _ = track(speed)
    lost = np.flatnonzero(e > R)
    print(f"speed {speed}: largest motion {step:4.1f} px/frame; error mean {e.mean():5.2f} "
          f"max {e.max():6.2f} px; mean iterations {n.mean():.1f}; "
          f"first frame with error > R: {lost[0] + 1 if len(lost) else 'none'}; "
          f"final Bhattacharyya coefficient {rho:.3f}")
e, _, _, c = track(7.0)
print("speed 7.0, error per frame (px):", np.round(e, 1).tolist())
print("speed 7.0, final estimate", np.round(c, 1), "; distractor centre (39.5, 112.0)")

speed 1.0: largest motion  3.0 px/frame; error mean  0.41 max   0.67 px; mean iterations 2.6; first frame with error > R: none; final Bhattacharyya coefficient 0.986
speed 3.0: largest motion  8.9 px/frame; error mean  0.80 max   4.55 px; mean iterations 4.3; first frame with error > R: none; final Bhattacharyya coefficient 0.989
speed 5.0: largest motion 14.8 px/frame; error mean  0.69 max   5.92 px; mean iterations 5.0; first frame with error > R: none; final Bhattacharyya coefficient 0.988
speed 7.0: largest motion 20.9 px/frame; error mean 43.43 max  87.36 px; mean iterations 2.9; first frame with error > R: 13; final Bhattacharyya coefficient 0.766
speed 9.0: largest motion 26.8 px/frame; error mean 52.27 max 102.44 px; mean iterations 1.4; first frame with error > R: 3; final Bhattacharyya coefficient 0.000
speed 7.0, error per frame (px): [0.2, 0.2, 0.6, 0.3, 0.5, 0.2, 0.4, 1.7, 3.1, 1.2, 0.9, 4.3, 26.2, 43.2, 59.2, 73.7, 83.1, 87.4, 85.3, 78.2, 68.1, 58.4, 52.0, 52.4, 55.7, 59.4, 58.8, 53.8, 45.7, 37.5, 37.0, 46.9, 61.2, 74.9, 82.3, 86.4, 82.7, 72.9, 58.1]
speed 7.0, final estimate [ 39.5 110.6] ; distractor centre (39.5, 112.0)

Up to a motion of about 15 pixels per frame, more than the disk’s radius, the tracker stays within a few pixels of the true centre, with an average error under one pixel and a few iterations per frame; the number of iterations grows with the motion, as the authors observed. At 21 pixels per frame the window keeps up for 12 frames and then locks onto the distractor’s colours: the error jumps and the final estimate sits within 1.5 pixels of the distractor’s centre with a Bhattacharyya coefficient of 0.77, a confident wrong answer. At 27 pixels per frame the disk escapes the window by the third frame, no pixel in the window has a target colour, and the tracker stops where it is. Both failures are the ones Comaniciu, Ramesh and Meer (2003, §7) warn about: motion larger than the basin of attraction calls for a motion model or several starting points, and the colour model must be discriminative enough for the scene.

Properties and Behavior

• Nonparametric. The number and shape of clusters come from the data and the bandwidth, not from a model; the results are deterministic for given seeds.
• Robust to outliers. Samples outside the window have no influence; Comaniciu and Meer (2002, §2.6) show that the iteration computing a robust location M-estimate is the mean shift iteration, and that the estimate is the mode of the density estimate.
• Convergence. Monotone ascent of the density estimate for convex, decreasing profiles, finite convergence with the flat kernel, and in general a linear rate that can be very slow near bandwidths where modes merge (see Convergence).
• Edge preservation. In the joint domain, filtering and segmentation keep discontinuities sharp, as the second example shows.
• Tracking. Fast and tolerant of deformation and partial occlusion because a histogram ignores the arrangement of colours, but weak when the background has the target’s colours, and limited to motions within the window’s reach.

Limitations

• One global bandwidth. A single $h$ cannot suit both a sharp, dense mode and a broad, sparse one; see Mean Shift Clustering for the condition Comaniciu and Meer (2002, §3) state.
• Scaling of features. The method assumes a Euclidean metric; features in different units must be rescaled or given their own bandwidths, as in the joint domain.
• Cost. All-pairs distances make naive mean shift quadratic in the number of samples per iteration, and convergence with smooth kernels can take many iterations.
• High dimensions. Kernel density estimates become unreliable as the dimension grows (see Complexity), which limits mean shift to low-dimensional feature spaces.
• Plateaus and saddle points. Iterations stop on flat regions and can, in rare cases, end at saddle points, so modes need post-processing; with coarse seeding, spurious stationary points can survive, as in the first example.
• Tracking assumptions. A colour histogram discards spatial layout and fails against similarly coloured backgrounds, and the basic tracker has no motion model.

Variants

• Blurring mean shift moves the data themselves at every iteration (Cheng, 1995), so the data set contracts toward its modes; it was the form Fukunaga and Hostetler (1975) proposed, as Carreira-Perpiñán (2007, §1) describes it.
• Gaussian mean shift as EM opens EM acceleration techniques, such as sparse EM, which updates the posterior weights of only a small set of nearby samples at most iterations, and switching to Newton’s method near a mode (Carreira-Perpiñán, 2007, §6).
• Adaptive bandwidths choose a bandwidth per sample or per location; Comaniciu and Meer (2002, §3.1) refer to their adaptive mean shift procedure built on a sample-point density estimator.
• Kernel-based tracking with scale adaptation and CamShift extend the procedure to tracking, as described above.
• Quick shift (Vedaldi & Soatto, 2008) also assigns each pixel to a mode of a density over colour and position and is used to compute superpixels.
• Other uses in vision apply mean shift to merge overlapping detections in sliding window detection, as Dalal’s thesis did in position and log-scale, and to cluster pixel embeddings into instances in instance segmentation.

History

Fukunaga and Hostetler (1975) studied the estimation of the gradient of a density and its uses in pattern recognition. By Carreira-Perpiñán’s account (2007, §1), they were perhaps the first to propose the idea and introduced the term mean shift, deriving it with the Epanechnikov kernel as gradient ascent on the log-density with a variable step size, without a convergence proof, and moving every data point at each iteration. Comaniciu and Meer (2002, §1) describe the procedure as largely forgotten until Cheng (1995) revived it; Cheng named the blurring process, discussed its convergence, and proposed iterating without moving the data. Comaniciu and Meer gave a convergence proof for kernels with convex profiles, whose step from the density values to the points was later found to be incomplete (see Convergence), and made the procedure a general tool for vision: colour clustering, and discontinuity-preserving filtering and segmentation in the joint domain, whose code they released; with Ramesh they added kernel-based tracking (2000; 2003). Carreira-Perpiñán (2007) connected Gaussian mean shift to EM and analysed its speed of convergence.