DBScan算法

本文依据《数据挖掘概念与技术》第三版（韩家炜）中的DBScan算法来实现的，其伪代码如下（截图）：

因为实验的需要，本文用经纬度数据作为测试数据，该距离计算参考了博主的文章（https://blog.csdn.net/xiejm2333/article/details/73297004）；如需要直角坐标数据，则将dbEpsilon函数中的两点距离距离计算公式直接改为欧拉公式即可。

算法中的数据结构说明：

实验所用测试数据地址：链接: https://pan.baidu.com/s/1JF74dKSWDum83XMozkjlPA 密码: 24fw

聚类结果保存在 C_all中，并以cell的形式对各簇加以区分；同时，cluster.txt也保存了结果，以0 0 0 0 加以区分各簇

算法的主程序：

clc;
clear;
%设定半径,公里
epsilon=100;
%邻域密度阈值，发生次数
minPts=10;
%读取数据文件,生成点矩阵
fileID = fopen('D:\matlabFile\GTD\gtd.txt');
 DS=textscan(fileID,'%f %f');
 fclose(fileID);
 %将cell类型转换为矩阵类型,这里的数据二维坐标点，附加其类型属性
 %组合DS中的前两列,由cell转换为矩阵
 Dataset=cat(2,DS{1},DS{2});
 %删除空值，默认值为999
index=find(Dataset(:,1)==999);
Dataset(index,:)=[];
countD=size(Dataset,1);
 %零矩阵
 temp=zeros(countD,1);
%增加第三列，全部置位0，代表unvisited，如果为1，代表visited，如果为2，代表噪声
 Dataset=cat(2,Dataset,temp);
  %生成整数序列,与每一个点对应
  sequence=zeros(1,countD);
sequence=(1:1:countD);
 %增加第四列，给坐标添加序号
 Dataset=cat(2,Dataset,sequence');
 %该临时数据集包含所有标记为unvisited的点
 Dataset_temp=Dataset;
 %绘制初始点图
%  scatter(Dataset(:,1),Dataset(:,2),'filled');
%簇的集合
C_all=cell(1,1);
%所有符合条件的点集合，集中存储，方便查询
C=zeros(1,4);
flag=0;
while flag==0
    %判断Dataset是否遍历完毕
    if size(Dataset_temp,1)~=0
        %随机选择一个对象p,在剩余的Dataset中
        randC=randperm(size(Dataset_temp,1),1);
        p=Dataset_temp(randC,:);
        %将Dataset_temp中的这个点剔除
        Dataset_temp(randC,:)=[];
        %将p从Dataset中剔除，之后要重新加入，保证Dataset始终是完整的
        Dataset_1=Dataset;        
        %在原数据库中找到p,标记为visited
        Dataset(p(1,4),3)=1;
        %标记p为visited
        p(1,3)=1;
        Dataset_1(p(1,4),:)=[];         
        %寻找在半径范围之内的对象
        %存放临时对象
        N=dbEpsilon(p,Dataset_1,epsilon);
        %新建的簇，用来保存符合条件的p
        C_temp=p;
        %判断p邻域内对象数量是否满足minPts的要求
        if size(N,1)>=minPts
            %遍历临时对象集N,寻找新的对象
            flag_1=0;
            count=0;
            while flag_1==0
                %如计数器与N的行数相同，则结束循环
                if count==size(N,1)
                    flag_1=1;
                else
                    count=count+1;
                    for i=1:size(N,1)
                        %判断抽取出的对象是否是unvisited，否则标记为visited
                        if N(i,3)==0
                            t=N(i,:);
                            t(1,3)=1;
                            N(i,3)=1;
                            Dataset(t(1,4),3)=1;
                            %Dataset_temp中删除标记为visited的点，使得Dataset_temp里的数据始终为unvisited
                            ind1=find(Dataset_temp(:,4)==t(1,4));
                            Dataset_temp(ind1,:)=[];
                            %将抽取的对象从Dataset_2中剔除
                            Dataset_2=Dataset_1;
                            ind2=find(Dataset_2(:,4)==t(1,4));
                            Dataset_2(ind2,:)=[];
                            %寻找所抽取对象epsilon范围内的对象
                            N_temp=dbEpsilon(t,Dataset_2,epsilon);
                            %判断抽取对象邻域范围内的对象数量是否满足minPts条件
                            if size(N_temp,1)>=minPts;
                                %判断N_temp中点是否已经存在于N，如不存在，则加入
                                for j=1:size(N_temp,1)
                                    ind3=find(N(:,4)==N_temp(j,4));
                                    if isempty(ind3)
                                        N=cat(1,N,N_temp(j,:));
                                    end
                                end
                            end
                            %判断抽取对象是否已经是簇的成员
                            ind4=find(C(:,4)==t(1,4));
                            if isempty(ind4)
                                %将t添加到点集合
                                C=cat(1,C,t);
                                
                                %判断抽取对象是否已经是簇C_temp的成员
                                ind5=find(C_temp(:,4)==t(1,4));
                                if isempty(ind5)
                                    C_temp=cat(1,C_temp,t);
                                end
                            end
                        end
                    end                    
                end                
            end
            if size(C_temp,1)>minPts
                C_all{size(C_all,1)+1,1}=C_temp;
            end
        else
            %在原Dataset中找到p,标记为噪声
            ind6=find(Dataset(:,4)==p(1,4));
            Dataset(ind6,3)=2;            
        end
    else
        flag=1;
    end
end
C_all(1,:)=[];

FID=fopen('cluster.txt','wt');
for k=1:size(C_all,1)
    for i=1:size(C_all{k,1},1)
        for j=1:size(C_all{k,1},2)-1            
            fprintf(FID, '%f ', C_all{k,1}(i,j));
        end
         fprintf(FID,'\n');    
    end  
    fprintf(FID, '%f %f %f %f ', [0 0 0 0]);
    fprintf(FID,'\n');
end
fclose(FID);

dbEpsilon函数实现如下：

function result=dbEpsilon(p_t,Dataset,epsilon)
    N=zeros(1,4);
    %地球半径，单位：公里
    R=6378.137;
    %角度转换成弧度
    p=[rad(p_t(1,1)),rad(p_t(1,2))];
    for i=1:size(Dataset,1)
        %计算两个对象的欧式距离
        q_temp=Dataset(i,:);
        q=[rad(q_temp(1,1)),rad(q_temp(1,2))];
        %两点纬度之差
        lat=p(1,1)-q(1,1);
        %两点经度之差
        lon=p(1,2)-q(1,2);
        %计算两点之间的距离
        dist=2*asin(sqrt((sin(lat/2.0))^2+cos(p(1,1))*cos(q(1,1))*(sin(lon/2.0))^2));
        %弧长乘以地球半径
        dist=dist*R;

        if dist<=epsilon
            N=cat(1,N,q_temp);
        end
    end
    %删除第一行空值
    N(1,:)=[];
    result=N;
end

function result=rad(value)
   result=value * pi/180.0; 
end