DBSCAN的理解和matlab实现

DBSCAN是基于密度的聚类算法，以下总结一下编写matlab时遇到的一些问题。

1、算法的基本流程

步骤1 ： 首先初始化变量，主要包括原始数据变量（此处为一个二维矩阵，包括x，y坐标，共1500个采样点），由randmperm生成的随机标签向量（一个一维的列向量），这个向量主要是用来随机挑选数据中的一个点开始分类。初始化数据的分类代号向量，这是一个一维列向量，他的值表征了每一个采样点是被分为哪一类。计算距离矩阵，这是一个二维的矩阵，主要是计算个点之间的距离。

步骤2 ：生成用于分类的数据，参考程序中的生成方法

步骤3：开始按照randmperm的随机标签来遍历数据集中的数据，若当前遍历到的数据点已经被分类完成（即其分类代号向量中的值不再是0，而被赋予了其他值），那么跳过已完成的分类点，如果当前点没有被分类过，那么进入下一个步骤。

步骤4：按照DBSCAN的想法，我们首先要搜寻半径为rmax领域内所有的点，搜索完成后得到领域内数据点的个数，以及他们在数据矩阵中的位置。之后判定领域内数据点的个数是否满足最小个数要求，如果不满足，那么这是个噪声点，直接将其分类代号赋值为-1（表示噪声点）；如果他满足领域内最小点数要求，那么进入下一步

步骤5：这一步关键要想清楚一件事前，假设我们首先初始化一个队列向量，他用来存储第4步中搜索得到的领域中所有点的标号，按照算法，下一步需要以领域中的每一个点在作为中心点去搜寻他们各自领域中的点，判定是否满足条件，这里有点别捏，我们详细的说明

假设初始化队列Qurey = [] ；%给他赋空值

并且 Qurey = [Qurey epsilonspace];%把第四步中遍历得到领域内的点首先给到队列之中，其中

                                                        %epsilonspace包含的是第四步中搜寻得到的领域内的所有点

//由于Qurey中的点还要在继续遍历，所以他的长度是要根据数据来自动增加的，每遍历一个Qurey中的点的领域，那些符合条件的点在按顺序补齐在Qurey之后，然后知道遍历完成所有的点，即结束，所以按照这个思路，实现的方法如下：

Qurey = [Qurey epsilonspace];   %把核心点在领域范围内的点都取出来
                      
countervalid = length(Qurey);
kcount = 1;

%开始对中心点内的所有点做遍历，由于此处需要不断地添加Qurey的个数，但是不同核心点中领 %域内可能存在相同的点，就是C点即在A的领域内，又在B的领域内，所以为了避免这种情况，需 %要删除在遍历过程中相同的点，
while(kcount <= countervalid) %这么写的目的就是为了方便删除相同的点

            Q = Qurey(kcount);

           %这是个领域搜寻函数，他的结果是输出领域内点的个数以及符合点的位置
            [tempnumber tempspace] = Searchepsilon(Q , rmax , clustter);

            %这个是个bug修复，因为距离矩阵计算时自己和自己的距离也算了是0
            tempspace(find(tempspace == Q)) = [];
            if tempnumber > pmin%pmin是聚类的最小要求点数
                difftemp = setdiff(tempspace , Qurey);
                Qurey = [Qurey difftemp];        %这个是把不同的点在重新添加到队列中
                clustter(Qurey) = clusttercount; %clusttercount表示了类的标号
            else
               clustter(Q) = clusttercount;   %表征这是边界点 ,任然归为当前类，此处不做特殊处理
            end
            kcount = kcount + 1;
            countervalid = length(Qurey);
        end
        clusttercount = clusttercount + 1;
    end 

步骤6 在完成步骤5后，相当于是一个点他完成了与他所有密度想通点的分类，我们可以跳转到第三步开始选择其他的点开始遍历。

下面是matlab代码

clearvars; close all;clc;
dbstop if error;
model_class = 3;
dim = 2;
% 期望值
m = [0, 0;
    2, 2;
    -2, -2];
% 协方差阵
s(:, :, 1) = [0.1, 0;
              0, 0.1];
s(:, :, 2) = [0.5, 0.3;
              0.3, 0.5];
s(:, :, 3) = [1, -0.5;
              -0.5, 1];
          
num = [500, 500, 500];
data = generate_data_GMM(dim, model_class, m, s, num);
data = data(:, (1:2));
figure;
scatter(data(:,1),data(:,2));
xlim([-5 5]);
ylim([-5 5]);

%% DBSCAN
global distanmatrix;
datasize = size(data,1);
index = randperm(datasize);
clustter = zeros(datasize , 1);     %数据样本点的分类信息 0为初始值，-1为噪声点值
distanmatrix = squareform(pdist(data));

rmax = 0.3 ; 
pmin = 12;
clusttercount = 1;

for i = 1 : 1 : datasize    %开始遍历数据库中的数据，是以随机序列的顺序来完成遍历的
    
    currentindex = index(i);
    currentpoint = data(currentindex,:);    %获取当前的数据
    
    if clustter(currentindex) ~= 0
       continue; 
    end
    Qurey = [];
    [epsilonnumber epsilonspace] = Searchepsilon(currentindex , rmax , clustter);  % 获取当前点的epsilon领域内的所有点,以及他们对应的标号
    
    if epsilonnumber < pmin
        clustter(currentindex) = -1;        %如果不满足领域个数条件，判定该点为噪声点    
    else
        clustter(currentindex) = clusttercount;
        epsilonspace(find(epsilonspace == currentindex)) = [];
        Qurey = [Qurey epsilonspace];   %把核心点在领域范围内的点都取出来
                      
        countervalid = length(Qurey);
        kcount = 1;
        cnt = 1;
        while(kcount <= countervalid)   %开始对中心点内的所有点做遍历，由于此处需要不断地添加Qurey的个数，所以查询完成一个就删除一个
            Q = Qurey(kcount);
            [tempnumber tempspace] = Searchepsilon(Q , rmax , clustter); 
            tempspace(find(tempspace == Q)) = [];
            if tempnumber > pmin
                difftemp = setdiff(tempspace , Qurey);
                if(~isempty(difftemp))                   
                    cnt = cnt + 1;
                end
                Qurey = [Qurey difftemp];
                clustter(Qurey) = clusttercount;
            else
               clustter(Q) = clusttercount;   %表征这是边界点 ,任然归为当前类，此处不做特殊处理
            end
            kcount = kcount + 1;
            countervalid = length(Qurey);
        end
        clusttercount = clusttercount + 1;
    end  
end

for i = 1 : 1 : clusttercount-1
   
    scatter(data(find(clustter == i),1),data(find(clustter == i),2));hold on;
end
scatter(data(find(clustter == -1),1),data(find(clustter == -1),2) , '*');hold on;
xlim([-5 5]);
ylim([-5 5]);
grid on;

% 领域搜寻函数，以currentpoint为中心点，搜寻其领域为rmax范围内所有的点的个数以及其领域内的所有点；
function [epsilonnumber epsilonspace] = Searchepsilon(currentindex , rmax , clustter)
    global distanmatrix;    %全局变量
    valiedindex = find(distanmatrix(currentindex , :) <= rmax);
%     deletindex = find(clustter ~= 0);
%     for i = 1 : length(deletindex)
%        
%         valiedindex(find(valiedindex == deletindex(i))) = [];
%     end
    

    epsilonnumber = length(valiedindex);
    epsilonspace = valiedindex;
end



%%
function data = generate_data_GMM(dim, model_class, m, s, num)
    %   生成多组多维高斯分布数据
    %   dim为高斯分布的维数
    %   model_class为生成的数据组数
    %   m为高斯分布的期望值，大小为model_class*dim
    %   s为高斯分布的协方差阵，大小为dim*dim*model_class
    %   num为各组高斯分布的数据量，其大小为model_class*1
    %   返回值data为生成的数据，其大小为num*(dim+1)
    %   前dim列为高斯分布数据，第(dim+1)列为组的类别编号
    %   该函数用于生成多组高斯分布数据，可为聚类算法提供数据
    data = [];
    for i = 1 : model_class
        data1 = ones(num(i), dim + 1);
        data1(:, (1 : dim)) = mvnrnd(m(i, :), s(:, :, i), num(i));
        for j = 1 : num(i)
            data1(j, dim + 1) = i;
        end
        data = [data; data1];
    end
end

 左边是原始数据，右边是分类完成的数据，“ * ”号点表示是噪声点，原点表示的是分类的点