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的控制可以参阅下图