基于密度的空间的数据聚类方法DBSCAN(Density-Based Spatial Clustering of Applications with Noise)matlab实现

 DBSCAN（Density-Based Spatial Clustering of Applications with Noise，具有噪声的基于密度的聚类方法）是一种基于密度的空间聚类算法。该算法将具有足够密度的区域划分为簇，并在具有噪声的空间数据库中发现任意形状的簇，它将簇定义为密度相连的点的最大集合。

 同一类别的样本，他们之间的紧密相连的，也就是说，在该类别任意样本周围不远处一定有同类别的样本存在。通过将紧密相连的样本划为一类，这样就得到了一个聚类类别。通过将所有各组紧密相连的样本划为各个不同的类别，则我们就得到了最终的所有聚类类别结果。

1. 代码解析

%% 函数功能注释
% Function: [class,type]=dbscan(X,k,Eps)
% -------------------------------------------------------------------------
% 目标: 
% 使用DBSCAN(Density-Based Scan Algorithm with Noise)聚类数据
% -------------------------------------------------------------------------
% 输入: 
% X - 数据集 (m,n); m-行，对象的个数, n-列，属性个数
% k - 对象周围邻居节点个数，一个聚类簇中最小的节点数
% Eps - ε半径，邻接点半径, 如果不清楚可以置为空
% -------------------------------------------------------------------------
% 输出:
% class - 向量，聚类簇
% type - 向量，第i个对象是什么类型
% (core: 1, border: 0, outlier: -1)
% 核心点：1 边界点：0 噪声点：-1
% -------------------------------------------------------------------------

%% 函数体
function [class,type]=dbscan(X,k,Eps)

[m,n]=size(X);

% 如果输入参数小于3个或输入半径为空
if nargin<3 || isempty(Eps)
   [Eps]=epsilon(X,k);
end


X=[[1:m]' X]; % 将X加上行序号1,2,3,4,5...
[m,n]=size(X);
type=zeros(1,m); % 初始化类型为0，均为边界值
no=1;
touched=zeros(m,1); % 记录节点有没有被计算过；0：没有计算过；1：计算过；


for i=1:m
    if touched(i)==0;
       ob=X(i,:);
       D=dist(ob(2:n),X(:,2:n)); % 取每一行的第二列到最后一列，即求当前对象与所有对象的距离
       ind=find(D<=Eps); % 查找距离小于给定半径对象的下标,如果有5个小于Eps的，则ind的长度为5
    
       % 边界点，即当前节点ε邻域内的节点，算上当前节点，所以小于k+1
       if length(ind)>1 && length(ind)<k+1       
          type(i)=0;
          class(i)=0;
       end
       
       % 噪声点，当前节点ε邻域内只有它本身
       if length(ind)==1
          type(i)=-1;             
          class(i)=-1;            
          touched(i)=1;
       end

       % 核心点，ε邻域内节点数大于k+1，算上当前节点，所以大于等于k+1
       if length(ind)>=k+1; 
          type(i)=1;
          class(ind)=ones(length(ind),1)*max(no); % 以核心点为中心的聚类簇，class里面放的是核心点的直接密度可达点的下标

          while ~isempty(ind)
                ob=X(ind(1),:); % 每次取第一行元素
                touched(ind(1))=1; % 被计算标记置为1
                ind(1)=[]; % 清除计算过的节点序号
                D=dist(ob(2:n),X(:,2:n)); % 取每一行的第二列到最后一列，即求当前对象与所有对象的距离
                i1=find(D<=Eps); % 查找距离小于给定半径对象的下标,如果有5个小于Eps的，则i1的长度为5
                
                % 边界点，即当前节点ε邻域内的节点，算上当前节点，所以小于k+1
                if length(i1)>1
                   class(i1)=no;
                   if length(i1)>=k+1;
                      type(ob(1))=1; % 核心点
                   else
                      type(ob(1))=0; % 边界点
                   end

                   for i=1:length(i1)
                       if touched(i1(i))==0
                          touched(i1(i))=1;
                          ind=[ind i1(i)];   
                          class(i1(i))=no;
                       end                    
                   end
                end
          end
          no=no+1; 
       end
   end
end

i1=find(class==0);
class(i1)=-1;
type(i1)=-1;

%% 函数功能注释
% function: [D]=dist(i,x)
%
% 目标: 
% 计算当前节点到所有节点的距离，返回一个距离向量D        
%                                                                    
% 输入: 
% i - 一个对象，1*n的向量
% x - 数据集 (m,n); m-行，对象的个数, n-列，属性个数          
%                                                                 
% 输出: 
% D - 欧氏距离，m*1
% -------------------------------------------------------------------------

%% 函数体
function [D]=dist(i,x)
    [m,n]=size(x);
    D=sqrt(sum((((ones(m,1)*i)-x).^2)'));

    if n==1
       D=abs((ones(m,1)*i-x))';
    end
end

%% 函数功能注释
% Function: [Eps]=epsilon(x,k)
%
% 目标: 
% 计算邻域半径
%
% 输入: 
% X - 数据集 (m,n); m-行，对象的个数, n-列，属性个数
% k - 对象周围邻居节点个数，一个聚类簇中最小的节点数
% -------------------------------------------------------------------------

%% 函数体
function [Eps]=epsilon(x,k)

    [m,n]=size(x);

    Eps=((prod(max(x)-min(x))*k*gamma(.5*n+1))/(m*sqrt(pi.^n))).^(1/n);
end

%% 函数功能注释
% Function: rundbscan(X,k,Eps)
% -------------------------------------------------------------------------
% 目标: 
% 调用dbscan函数
% -------------------------------------------------------------------------
% 输入: 
% X - 数据集 (m,n); m-行，对象的个数, n-列，属性个数
% k - 对象周围邻居节点个数，一个聚类簇中最小的节点数
% Eps - ε半径，邻接点半径, 如果不清楚可以置为空
% -------------------------------------------------------------------------

%% 函数体
function rundbscan(X,k,Eps)
    [c,t] = dbscan(X,k,Eps);
    [m,n] = size(X);
    figure
    axis equal
    hold on;
    for i = 1:m
        switch t(i)
            case 1
                % 核心点为红色 并以核心节点为圆心，ε为半径画圆
                plot(X(i,1),X(i,2),'ro');
                drawCircle(X(i,1),X(i,2),0.2);
            case 0
                plot(X(i,1),X(i,2),'go');
            case -1
                plot(X(i,1),X(i,2),'bo');
        end
    end
end

测试数据：X的值

    0.5604    0.8121
    0.5978    0.1031
    0.1921    0.5163
    0.2263    0.1556
    0.1540    0.2788
    0.3778    0.8826
    0.1120    0.6316
    0.5113    0.9637
    0.6385    0.5997
    0.1222    0.0095
    0.3442    0.4156
    0.0703    0.0224
    0.5258    0.3140
    0.6503    0.4816
    0.6155    0.5755
    0.8139    0.2387
    0.2356    0.5583
    0.6182    0.3949
    0.3028    0.5066
    0.0888    0.1960
    0.3053    0.1568
    0.0714    0.7545
    0.9586    0.3300
    0.7272    0.6442
    0.1015    0.5655
    0.7863    0.8995
    0.3522    0.8807
    0.9114    0.4431
    0.4955    0.2047
    0.8727    0.6322
    0.2205    0.6368
    0.3672    0.6924
    0.3668    0.9075
    0.0971    0.7826
    0.3310    0.6253
    0.6597    0.6010
    0.4665    0.2128
    0.8792    0.6514
    0.2697    0.8145
    0.3548    0.9223
    0.3204    0.9167
    0.5517    0.5611
    0.1882    0.2598
    0.1732    0.3338
    0.9395    0.9978
    0.3766    0.5256
    0.5871    0.8159
    0.0552    0.5970
    0.1409    0.5069
    0.0964    0.3485

运行结果：

设置聚类半径： Eps = 0.19；聚类密度：k = 8；

红色为核心点，绿色为边界点，蓝色为噪声点；

2. 参考资料

关于DBSCAN的原理可以参考如下几篇文章，写的都挺好，看完这几篇文章再看我的代码。

https://www.cnblogs.com/pinard/p/6208966.html

https://www.cnblogs.com/chaosimple/archive/2013/07/01/3164775.html

http://blog.csdn.net/zhouxianen1987/article/details/68945844