DBSCAN聚类分析在图像分割的应用

工具：matlab

k均值和dbscan我就不说了，网上很容易就搜到了。
我这个是基于k均值的dbscan，是当时老师给我发的一篇论文，我看了之后把东西稍微改了一下拿来用的，因为原来的那个论文的算法里面有个动态参数他没说是怎么得到的，然后我就把它变了，直接设成常数了，因为涉及到的因变量的值是两个方向决定的，一个设成常数，通过另一个来体现就可以了。
基本流程就是先把数据读进去，然后改成灰度值数据，再进行聚类，完了再给他染色。
步骤1：先将原始数据进行K均值聚类处理，得到含有 个标记的矩阵；
步骤2：对第 个类求出类内最大距离和类内最小距离；
步骤3：再将最大距离减去最小距离分为 个组，构建一个统计量，每段距离下面对应这个距离包含了样本对数（一对样本产生一个距离）；
步骤4：找出包含样本对数最多的一组距离，取他们的距离的中值作为该类的距离参数 （比如最大距离20，最小距离5， 是5，就用20减去5得到15,再将距离15分成5段，每一段占3，第一段就是6到8，第二段9到11，第三段12到14，第四段15到17，第五段18到20。假如第四段距离包含了最多的样本对数，就将（15+17）/2得到的值16作为该类的距离 ）；
步骤5：计算类中每个元素的r邻域包含的元素个数，取最小的作为这个类的数量参数 。所有类的数量参数确定之后，选取最大的作为所有类的数量参数 ；
步骤6：现在一共有 个距离参数 和一个数量参数 。重新对数据进行DBSCAN聚类。选择一个没考虑过的元素，将集合 中所有与之密度相连 的元素找出，归为一类，如果没有，则将其标为噪声点；
步骤7：重复步骤1直到 中所有的元素都被考虑过。
本来想做RGB数据的聚类的，三维的，灰度值数据是一维的，换成三维后需要的内存太大了！直接运行不下去出错了，然后就放弃了，要是不想放弃就改个距离，不用欧氏距离。
结果：





足球选手图像是我同学给我发的，哈哈。那天聊天说了一下我的进度，他就拿来让我试试，试试就逝世。
Im_sage_nospace:（主函数?，不知道怎么叫，就是不被调用的那种）
先看一下直方图，自己选个k在下面赋值，然后再运行。

I=imread('C:\Users\Administrator\Desktop\图像分割的图\the.jpg');
I1=rgb2gray(I);
subplot(2,2,1);imshow(I);title('原始图像')
subplot(2,2,2);imshow(I1);title('灰度图')
subplot(2,2,3);imhist(I1);
 
K=3;
epsilon=5;
IDX=DBSCAN_kmeans(I1,K,epsilon);
P=unique(IDX);
subplot(2,2,4);imshow(label2rgb(IDX));title('分割图')

DBSCAN_kmeans:

function IDX=DBSCAN_kmeans(GRAY,K,epsilon)
    [m1,m2]=size(GRAY);
    X=reshape(GRAY,m1*m2,1);
    X=double(X);
    b=kmeans(X,K);
    MinPts=zeros(K,1);
    
    for i =1:K
        G=X;
        G((b(:)~=i))=[];
        G=unique(G);
        G(G(:)==0)=[];
        y=zeros(numel(G),1);
        for q=1:numel(G)
         y(q)=numel(RegionQuery(G(q)));
        end
        MinPts(i)=sum(y)/numel(y); 
    end
       
    C=0;
    n=size(X,1);
    IDX=zeros(n,1);
    visited=false(n,1);   
    for i=1:n
        if ~visited(i)
            visited(X==X(i))=true;
            l=b(i);
            Neighbors=RegionQuery(X(i));                     
            if numel(Neighbors)>=MinPts(l) 
                B=zeros(2*n,1);
                B(1:numel(Neighbors),1)=Neighbors;
                C=C+1;
                ExpandCluster(i,C);
            end            
        end
    end
    IDX=reshape(IDX,m1,m2);
    
    function ExpandCluster(i,C)        
        IDX((X<=X(i)+epsilon)&(X>=X(i)-epsilon))=C;       
        visited(X==X(i))=true;
        k = 1;
        while true            
            k=k+1;          
            j = B(k);    
             if j==0
                break;
            end
            if ~visited(j)
                visited(j)=true;
                Neighbors2=RegionQuery(X(j));                
                visited(X==X(j))=true;
                l=b(j);
                if numel(Neighbors2)>=MinPts(l)
                    IDX((X<=X(j)+2)&(X>=X(j)-2))=C;
                    t=sum(sum(B~=0));                                           
                    B(t+1:t+numel(Neighbors2),1)=Neighbors2;
                    B1=unique(B);
                    B1(1)=[];
                    B=zeros(2*n,1);
                    B(1:numel(B1),1)=B1;
                end
            end            
                 
        end
    end
    
    function Neighbors=RegionQuery(i)
        D=abs(X-i);
        Neighbors=find(D<=epsilon);
    end
end

下面这个是灰度变成RGB的，我是从csdn上找到一个，然后也给里面改了改，改成合适我的，就拿来用了

function R=gray2rgb(img1)
% img1 - Source Image  (gray image)   
% img2 - Selected color image for coloring the gray image. 
clc;
warning off;
imt=img1;
ims=imread('C:\Users\Administrator\Desktop\图像分割的图\qise.png');
[sx, sy, sz]=size(imt);
[tx, ty ,tz]=size(ims);
if sz~=1
    imt=rgb2gray(imt);
end
if tz~=3
    disp ('img2 must be a color image (not indexed)');
else
    imt(:,:,2)=imt(:,:,1);
    imt(:,:,3)=imt(:,:,1);

% Converting to ycbcr color space
    nspace1=rgb2ycbcr(ims);
    nspace2= rgb2ycbcr(imt);
    ms=double(nspace1(:,:,1));
    mt=double(nspace2(:,:,1));
    m1=max(max(ms));
    m2=min(min(ms));
    m3=max(max(mt));
    m4=min(min(mt));
    d1=m1-m2;
    d2=m3-m4;
% Normalization
    dx1=ms;
    dx2=mt;
    dx1=(dx1*255)/(255-d1);
    dx2=(dx2*255)/(255-d2);
    [mx,my,mz]=size(dx2);
%Luminance Comparison
    disp('Please wait..................');
    for i=1:mx
        for j=1:my
             iy=dx2(i,j);
             tmp=abs(dx1-iy);
             ck=min(min(tmp));
             [r,c] = find(tmp==ck);
             ck=isempty(r);
             if (ck~=1)            
                 nimage(i,j,2)=nspace1(r(1),c(1),2);
                 nimage(i,j,3)=nspace1(r(1),c(1),3);
                 nimage(i,j,1)=nspace2(i,j,1);           
            end
         end
     end
    rslt=ycbcr2rgb(nimage);
    R=uint8(rslt);
end
end