Meanshift，聚類演算法
http://www.cnblogs.com/liqizhou/archive/2012/05/12/2497220.html
記得剛讀研究生的時候，學習的第一個演算法就是meanshift演算法，所以一直記憶猶新，今天和大家分享一下Meanshift演算法，如有錯誤，請線上交流。
Mean Shift演算法,一般是指一個迭代的步驟,即先算出當前點的偏移均值,移動該點到其偏移均值,然後以此為新的起始點,繼續移動,直到滿足一定的條件結束.
1. Meanshift推導
給定d維空間Rd的n個樣本點 ,i=1,…,n,在空間中任選一點x，那麼Mean Shift向量的基本形式定義為:
Sk是一個半徑為h的高維球區域,滿足以下關係的y點的集合,
k表示在這n個樣本點xi中,有k個點落入Sk區域中.
以上是官方的說法，即書上的定義，我的理解就是，在d維空間中，任選一個點，然後以這個點為圓心，h為半徑做一個高維球，因為有d維，d可能大於2，所以是高維球。落在這個球內的所有點和圓心都會產生一個向量，向量是以圓心為起點落在球內的點位終點。然後把這些向量都相加。相加的結果就是Meanshift向量。
如圖所以。其中黃色箭頭就是Mh（meanshift向量）。
再以meanshift向量的終點為圓心，再做一個高維的球。如下圖所以，重複以上步驟，就可得到一個meanshift向量。如此重複下去，meanshift演算法可以收斂到概率密度最大得地方。也就是最稠密的地方。
最終的結果如下：
Meanshift推導：
把基本的meanshift向量加入核函式，核函式的性質在這篇部落格介紹：http://www.cnblogs.com/liqizhou/archive/2012/05/11/2495788.html
那麼，meanshift演算法變形為
(1)
解釋一下K()核函式，h為半徑，Ck,d/nhd 為單位密度，要使得上式f得到最大，最容易想到的就是對上式進行求導，的確meanshift就是對上式進行求導.
(2)
令：
K(x)叫做g(x)的影子核，名字聽上去聽深奧的，也就是求導的負方向，那麼上式可以表示
對於上式，如果才用高斯核，那麼，第一項就等於fh,k
第二項就相當於一個meanshift向量的式子：
那麼（2）就可以表示為
下圖分析
的構成，如圖所以，可以很清晰的表達其構成。
要使得=0，當且僅當
=0，可以得出新的圓心座標：
（3）
上面介紹了meanshift的流程，但是比較散，下面具體給出它的演算法流程。
- 選擇空間中x為圓心，以h為半徑為半徑，做一個高維球，落在所有球內的所有點xi
- 計算，如果<ε(人工設定)，推出程式。如果>ε, 則利用（3）計算x，返回1.
2.meanshift在影象上的聚類：
真正大牛的人就能創造演算法，例如像meanshift，em這個樣的演算法，這樣的創新才能推動整個學科的發展。還有的人就是把演算法運用的實際的運用中，推動整個工業進步，也就是技術的進步。下面介紹meashift演算法怎樣運用到影象上的聚類核跟蹤。
一般一個影象就是個矩陣，畫素點均勻的分佈在影象上，就沒有點的稠密性。所以怎樣來定義點的概率密度，這才是最關鍵的。
如果我們就算點x的概率密度，採用的方法如下：以x為圓心，以h為半徑。落在球內的點位xi 定義二個模式規則。
（1）x畫素點的顏色與xi畫素點顏色越相近，我們定義概率密度越高。
（2）離x的位置越近的畫素點xi，定義概率密度越高。
所以定義總的概率密度，是二個規則概率密度乘積的結果，可以（4）表示
（4）
其中：
代表空間位置的資訊，離遠點越近，其值就越大，
表示顏色資訊，顏色越相似，其值越大。如圖左上角圖片，按照（4）計算的概率密度如圖右上。利用meanshift對其聚類，可得到左下角的圖。
|
|
|
|
|
|
----------------------------------------------------------------------------------------------------------------------
meanshift演算法
http://blog.sina.com.cn/s/blog_6fd15d5f01016agj.html
Mean Shift演算法,一般是指一個迭代的步驟,即先算出當前點的偏移均值,移動該點到其偏移均值,然後以此為新的起始點,繼續移動,直到滿足一定的條件結束.
1. Meanshift推導
給定d維空間Rd的n個樣本點 ,i=1,…,n,在空間中任選一點x，那麼Mean Shift向量的基本形式定義為:
Sk是一個半徑為h的高維球區域,滿足以下關係的y點的集合,
k表示在這n個樣本點xi中,有k個點落入Sk區域中.
以上是官方的說法，即書上的定義，我的理解就是，在d維空間中，任選一個點，然後以這個點為圓心，h為半徑做一個高維球，因為有d維，d可能大於2，所以是高維球。落在這個球內的所有點和圓心都會產生一個向量，向量是以圓心為起點落在球內的點位終點。然後把這些向量都相加。相加的結果就是Meanshift向量。
如圖所以。其中黃色箭頭就是Mh（meanshift向量）。
再以meanshift向量的終點為圓心，再做一個高維的球。如下圖所以，重複以上步驟，就可得到一個meanshift向量。如此重複下去，meanshift演算法可以收斂到概率密度最大得地方。也就是最稠密的地方。
最終的結果如下：
Meanshift推導：
把基本的meanshift向量加入核函式，核函式的性質在這篇部落格介紹：http://www.cnblogs.com/liqizhou/archive/2012/05/11/2495788.html
那麼，meanshift演算法變形為
(1)
解釋一下K()核函式，h為半徑，Ck,d/nhd 為單位密度，要使得上式f得到最大，最容易想到的就是對上式進行求導，的確meanshift就是對上式進行求導.
(2)
令：
K(x)叫做g(x)的影子核，名字聽上去聽深奧的，也就是求導的負方向，那麼上式可以表示
對於上式，如果才用高斯核，那麼，第一項就等於fh,k
第二項就相當於一個meanshift向量的式子：
那麼（2）就可以表示為
下圖分析的構成，如圖所以，可以很清晰的表達其構成。
要使得=0，當且僅當=0，可以得出新的圓心座標：
（3）
上面介紹了meanshift的流程，但是比較散，下面具體給出它的演算法流程。
- 選擇空間中x為圓心，以h為半徑為半徑，做一個高維球，落在所有球內的所有點xi
- 計算，如果<ε(人工設定)，推出程式。如果>ε, 則利用（3）計算x，返回1.
2.meanshift在影象上的聚類：
真正大牛的人就能創造演算法，例如像meanshift，em這個樣的演算法，這樣的創新才能推動整個學科的發展。還有的人就是把演算法運用的實際的運用中，推動整個工業進步，也就是技術的進步。下面介紹meashift演算法怎樣運用到影象上的聚類核跟蹤。
一般一個影象就是個矩陣，畫素點均勻的分佈在影象上，就沒有點的稠密性。所以怎樣來定義點的概率密度，這才是最關鍵的。
如果我們就算點x的概率密度，採用的方法如下：以x為圓心，以h為半徑。落在球內的點位xi 定義二個模式規則。
（1）x畫素點的顏色與xi畫素點顏色越相近，我們定義概率密度越高。
（2）離x的位置越近的畫素點xi，定義概率密度越高。
所以定義總的概率密度，是二個規則概率密度乘積的結果，可以（4）表示
（4）
其中：代表空間位置的資訊，離遠點越近，其值就越大，表示顏色資訊，顏色越相似，其值越大。如圖左上角圖片，按照（4）計算的概率密度如圖右上。利用meanshift對其聚類，可得到左下角的圖。
|
|
|
|
|
|
http://www.cnblogs.com/liqizhou/archive/2012/05/12/2497220.html
Matlab中meanshift演算法
mean-shift 的特點是把支撐空間和特徵空間在資料密度的框架下綜合了起來。對影象來講，支撐空間就是畫素點的座標，特徵空間就是對應畫素點的灰度或者RGB三分量。將這兩個空間綜合後，一個資料點就是一個5維的向量：[x,y,r,g,b]。
這在觀念上看似簡單，實質是一個飛躍，它是mean-shift方法的基點。
mean-shift方法很寶貴的一個特點就是在這樣迭代計算的框架下，求得的mean-shift向量必收斂於資料密度的區域性最大點。可以細看[ComaniciuMeer2002]的文章。
寫了點程式，可以對影象做簡單的mean-shift filtering，供參考：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DRGB, DSD, MSSD] = MScut(sMode, RGB_raw, hs, hf, m );
% designed for segmenting a colour image using mean-shift [ComaniciuMeer 2002]
% image must be color
% procedure in mean-shift
% 1. combine support space and feature space to make a mean-shift space
% based data description
% 2. for every mean-shift space data
% 3. do mean-shift filtering
% until convergence
% 4. end
% 5. find the converged mean-shift space data that you are interested in
% and label it
% 6. repeat the above steps
%
% a -- data in support space
% b -- data in feature space
% x -- data in mean-shift space
% f(.) -- data density function
% k(.) -- profile function (implicit)
% g(.) -- profile function (explicit)
% m -- mean shift vector
% hs -- bandwidth in support space
% hf -- bandwidth in feature space
% M -- threshold to make a distinct cluster
%% enter $hs$, $hf$, $m$ if necessary
if ~exist('hs')
hs = input('please enter spatial bandwidth (hs):n');
end
if ~exist('hf')
hf = input('please enter feature bandwidth (hf):n');
end
if ~exist('m')
m = input('please enter minimum cluster size (m):n');
end
switch upper(sMode)
case 'RGB'
RGB = double( RGB_raw );
case 'gray'
error('FCMcut must use colored image to do segmentation!')
end
sz = size(RGB);
mTCUT = Tcut( RGB(:,:,1) ); % trivial segmentation
%% project data into mean-shift space to make $MSSD$ (mean-shift space data)
mT = repmat([1:sz(1)]', 1, sz(2));
vX = mT(1:end)'; % row
mT = repmat([1:sz(2)], sz(1), 1);
vY = mT(1:end)'; % column
mT = RGB(:,:,1);
vR = mT(1:end)'; % red
mT = RGB(:,:,2);
vG = mT(1:end)'; % green
mT = RGB(:,:,3);
vB = mT(1:end)'; % blue
MSSD = [vX, vY, vR, vG, vB];
%% make $g$ - explicit profile function
disp('Using flat kernel: Epanechnikov kernel...')
g_s = ones(2*hs+1, 2); % 's' for support space
g_f = ones(2*hf+1, 3); % 'f' for feature space
%% main part $$
nIteration = 4;
nData = length(MSSD); % total number of data
DSD = MSSD*0; % 'DSD' for destination space data
for k = 1:nData
%
tMSSD = MSSD(k,:); % 't' for temp
for l = 1:nIteration
%
mT = abs( MSSD - repmat(tMSSD, nData, 1));
vT = logical( (mT(:,1)<=hs).*(mT(:,2)<=hs).*(mT(:,3)<=hf).*(mT(:,4)<=hf).*(mT(:,5)<=hf) );
v = MSSD(vT,:);
% update $tMSSD$
tMSSD = mean( v, 1 );
if nIteration == l
DSD(k,:) = tMSSD;
end
end
end
% show result
DRGB = RGB * 0;
DRGB(:,:,1) = reshape(DSD(:,3), sz(1), sz(2)); % red
DRGB(:,:,2) = reshape(DSD(:,4), sz(1), sz(2)); % red
DRGB(:,:,3) = reshape(DSD(:,5), sz(1), sz(2)); % red
figure, imshow(uint8(DRGB), [])
-------------------------------------------------------------------------------------------------------------------------
matlab練習程式（meanshift影象聚類）
http://www.cnblogs.com/tiandsp/archive/2012/11/20/2779601.html
關於這個meanshift，一來可以用來作為目標跟蹤，二來可以用來進行影象聚類。我這裡只實現了影象聚類，當然，是按自己的理解編寫的程式。至於目標跟蹤將來一定也是要實現的，因為我最初看這個演算法的原因就是想用他來跟蹤目標的。
meanshift的基本原理我就不介紹了，比起我的介紹，網上有不少牛人們比我解釋的好，最後我會列出我參考的文章。我這裡說一下我是怎麼理解meanshift影象聚類的。這裡的聚類也像過去的濾波一樣，需要一個模板矩陣，不過這個模板不是事先設定好的矩陣，而是在當前處理的畫素周圍提取一個r*r的矩陣，然後把這個矩陣化為一維向量，再對這個向量進行meanshift，最終迭代到的值再賦值給當前處理的畫素。所以可以這樣理解，把影象經過meanshift迭代到相同值的畫素聚為一類。
我這裡使用的是灰度影象，至於彩色影象，我看到一篇部落格上把rgb域轉換到luv域上再去做處理，這個我就不太清楚了，不過我看他的程式碼其中有一部分很像均值濾波。雖然我沒有和他用一樣的方法，不過他的程式碼也可以參考一下。傳送門在此。
下面是程式碼（這都是我自己的理解，不能保證都正確，不過至少可以為你的編碼提供一些思路）：
main.m
clear all;
close all;
clc;
r=2; %濾波半徑
img=imread('lena.jpg');
imshow(img);
img=double(img);
[m n]=size(img);
imgn=zeros(m+2*r+1,n+2*r+1);
imgn(r+1:m+r,r+1:n+r)=img;
imgn(1:r,r+1:n+r)=img(1:r,1:n);
imgn(1:m+r,n+r+1:n+2*r+1)=imgn(1:m+r,n:n+r);
imgn(m+r+1:m+2*r+1,r+1:n+2*r+1)=imgn(m:m+r,r+1:n+2*r+1);
imgn(1:m+2*r+1,1:r)=imgn(1:m+2*r+1,r+1:2*r);
imshow(mat2gray(imgn))
for i=1+r:m+r
for j=1+r:n+r
ser=imgn(i-r:i+r,j-r:j+r);
ser=reshape(ser,[1 (2*r+1)^2]); %將二維模板變為一維
imgn(i,j)=mean_shift(ser,2*r^2+2*r+1); %取模板最中間的那個值作為迭代初值
end
end
figure;
imgn=imgn(r+1:m+r,r+1:n+r);
imshow(mat2gray(imgn));
meanshift.m
function re= mean_shift( ser,p)
[m n]=size(ser);
tmp=double(ser);
pre_w=tmp(p);
point=p;
while 1
ser=tmp-pre_w;
for i=1:m*n
if i ~= point
ser(i)=ser(i)/(i-point); %i-point是距離，就是各種公式裡的h
end
end
ser=ser.^2;
K=(1/sqrt(2*pi))*exp(-0.5*ser); %傳說中的核函式
w=sum(tmp.*(K))/sum(K);
if abs(w-pre_w)<0.01
break;
end
pre_w=w;
end
% tmp1=abs(tmp-w);
% [i point]=min(tmp1);
re=w;
% if max(tmp)-w<0.01
% point=0;
% end
% point=w;
end
處理的效果：
原圖
半徑為2處理的效果
——————————下面是2013.5.30新增————————————
上一部分的meanshift影象聚類還需修改，下面實現最簡單的meanshift演算法，完全按照原理來。
最後的參考文獻都是很好的總結，不過這次我是參考的《影象處理、分析與機器視覺（第3版）》這本書。
下面是通常所見的迭代效果：
程式如下：
clear all; close all; clc;
%測試資料
mu=[0 0]; %均值
S=[30 0;0 35]; %協方差
data=mvnrnd(mu,S,300); %產生300個高斯分佈資料
plot(data(:,1),data(:,2),'o');
h=3; %核的大小
x=[data(1,1) data(1,2)]; %以第一個資料為迭代初值
pre_x=[0 0];
hold on
while norm(pre_x-x)>0.01;
pre_x=x;
plot(x(1),x(2),'r+');
u=0; %分子累加項
d=0; %分母累加項
for i=1:300
%最關鍵的兩步，均值位移公式實現
k=norm((x-data(i,:))/h).^2;
g=(1/sqrt(2*pi))*exp(-0.5*k);
u=data(i,:)*g+u;
d=g+d;
end
M=u/d; %迭代後的座標位置
x=M;
end
參考：
1.http://en.wikipedia.org/wiki/Mean-shift wiki百科，介紹的簡介明瞭。
2.http://www.cnblogs.com/liqizhou/archive/2012/05/12/2497220.html 非常詳細的理解。
3.http://emuch.net/bbs/viewthread.php?tid=4626864 小木蟲上一個同學的理解。
4.http://en.wikipedia.org/wiki/Kernel_(statistics) 介紹核函式的。
5.http://wenku.baidu.com/view/11b6a7de6f1aff00bed51eac.html 提出meanshift演算法的論文，雖然我沒怎麼看，不過想對演算法徹底理解的還是看這篇好。
------------------------------------------------------------------------------------------------------------
Meanshift影象平滑之opencv實現
http://www.cnblogs.com/easymind223/archive/2012/07/03/2574887.html
一句話一幅圖理解meanshift演算法：
對於集合中的每一個元素，對它執行下面的操作：把該元素移動到它鄰域中所有元素的特徵值的均值的位置，不斷重複直到收斂。
準確的說，不是真正移動元素，而是把該元素與它的收斂位置的元素標記為同一類。對於影象來說，所有元素程矩陣排列，特徵值便是畫素的灰度值。
Meanshift的這種思想可以應用於目標跟蹤、影象平滑、邊緣檢測、聚類等，是一種適應性很好的演算法，缺點是速度非常慢。
本文以影象平滑為例對其說明
從網上找程式碼不如自己動手寫。說明一下兩個引數的含義，hs和hr是核函式的視窗大小，hs是距離核函式，控制子視窗的大小，同時也影響計算速度。hr是顏色核函式，是顏色差值的閾值，maxiter是最大迭代次數。轉載請註明出處，謝謝。本文演算法只是用作實驗之用，沒有進行優化，計算時會有重複計算的地方，速度非常慢，且只支援3通道影象。
1 void MyTreasureBox::MeanShiftSmooth(const IplImage* src, IplImage* dst, int hs, int hr, int maxIter)
2 {
3 if(!src)return ;
4
5 IplImage* srcLUV = cvCreateImage( cvGetSize( src ), src->depth, src->nChannels );
6 IplImage* dstLUV = cvCreateImage( cvGetSize( src ), src->depth, src->nChannels );
7
8 cvCvtColor( src, srcLUV, CV_RGB2Luv);
9 cvCopy( srcLUV, dstLUV );
10
11 int widthstep = srcLUV->widthStep;
12 int channel = srcLUV->nChannels;
13
14 for( int y = 0; y<src->height; y++ )
15 {
16 for( int x = 0; x<src->width; x++ )
17 {
18 uchar L = (uchar)srcLUV->imageData[y *widthstep + x *channel];
19 uchar U = (uchar)srcLUV->imageData[y *widthstep + x *channel + 1];
20 uchar V = (uchar)srcLUV->imageData[y *widthstep + x *channel + 2];
21 int xx = x;
22 int yy = y;
23
24 int nIter = 0;
25 int count, sumL, sumu, sumv, sumx, sumy;
26
27 while(nIter < maxIter)
28 {
29 count = 0;
30 sumL = sumu = sumv = 0;
31 sumx = sumy = 0;
32
33 for( int m = y - hs; m <= y + hs; m++ )
34 {
35 for( int n = x - hs; n <= x + hs; n++ )
36 {
37 if(m >= 0 && m < src->height && n >= 0 && n < src->width)
38 {
39 uchar l = (uchar)srcLUV->imageData[m *widthstep + n *channel];
40 uchar u = (uchar)srcLUV->imageData[m *widthstep + n *channel + 1];
41 uchar v = (uchar)srcLUV->imageData[m *widthstep + n *channel + 2];
42
43 double dist = sqrt( (double)((L - l)^2 + (U - u)^2 + (V - v)^2) );
44 if( dist < hr )
45 {
46 count++;
47 sumL += l;
48 sumu += u;
49 sumv += v;
50 sumx += n;
51 sumy += m;
52 }
53 }
54 }
55 }
56 if(count == 0)break;
57 L = sumL / count;
58 U = sumu / count;
59 V = sumv / count;
60 xx = sumx / count;
61 yy = sumy / count;
62
63 nIter++;
64 }
65 dstLUV->imageData[y *widthstep + x *channel] = L;
66 dstLUV->imageData[y *widthstep + x *channel + 1] = U;
67 dstLUV->imageData[y *widthstep + x *channel + 2] = V;
68 }
69 }
70
71 cvCvtColor( dstLUV, dst, CV_Luv2RGB );
72 cvReleaseImage(&srcLUV);
73 cvReleaseImage(&dstLUV);
74 }
hs和hr的控制可以參閱下圖
