DBSCAN聚类算法的理解与应用

在前面的文章中，我们讲了KNN算法的原理与简单应用，KNN一种有监督学习的分类算法，也就是说该算法首先需要训练数据来进行学习之后才能对数据进行分类。在本文中我们讲到的DBSCAN聚类算法，也属于一种数据分类算法，只不过该算法不需要任何训练数据就能对数据进行分类，因此该算法属于无监督的数据分类算法。本文中我们首先讲一下该算法的原理，然后举一个例子来说明该算法的应用。

1. DBSCAN算法原理

首先介绍该算法的主要概念与参数：

(1) ε值：样本与样本之间的距离阈值，如果样本A与样本B的距离小于该阈值，则认为样本A在样本B的邻域内，同时样本B也在样本A的邻域内。

(2) minPts：每一个样本的邻域内样本数阈值，如果该样本邻域内的样本数大于等于该阈值，则认为该样本是核心点。

(3) 核心点：即邻域内的样本数大于等于minPts的样本。如下图所示，如果样本A的邻域内（以A为圆心的圆内）样本数达到minPts以上，则认为A为核心点。

(4) 样本距离：欧式距离与曼哈顿距离是两种很常见的衡量数据样本距离的指标，假设有样本A(a1,a2,...,an)和样本B(b1,b2,...,bn)，那么A与B的欧式距离为：

曼哈顿距离为：

(5) 样本的访问标记：一开始将所有样本的标记设置为-1，表示所有样本都没有被访问。算法执行过程中，会遍历一遍所有样本，经过遍历的样本则将其标记置1，表示该样本已经被访问过，不用再处理。

(6) 样本的类编号：设置一个初始类编号为-1，分类过程中，每新增一个类，类编号加1。每当一个样本被归类到某一个类之后，该样本的类编号则设置为当前新增的类编号。所以可以通过判断该样本的类编号是否为-1来判断其是否已经被归类。

DBSCAN算法的核心思想是，判断每一个样本是不是核心点，如果是核心点，则对该样本及其邻域内的样本进行分类处理，具体处理流程如下：

下面我们使用C++和opencv实现DBSCAN算法，对灰度图像进行分类。样本就是灰度图像的每一个像素点，可以把每一个像素点看成一个三维向量，由x坐标、y坐标、像素值I(x,y)组成[x, y, I(x,y)]。为了简化计算，我们使用曼哈顿距离来计算各样本之间的距离。

为了减少计算耗时，首先计算出所有样本之间的距离，并统计每个样本邻域内的样本数，以及每个样本邻域内所有样本的序号，方便后续的查找。下面上代码。

样本的类：

class point
{
  public:
    float x;    //x坐标
    float y;    //y坐标
    float z;    //像素值
    int cluster;    //簇的序号，为-1则表示不属于任何簇
    int pointType;   //1 噪点 2 边界点 3 核心点 
    int visited;     // 1 - 已访问   0 - 未访问
    int pts;        //该点周围邻域内的点数
    point()     //析构函数初始化
    {
      cluster = -1;
      pointType = 1;
      visited = 0;
      pts = 0;
    }
};

邻域内点的类：

class neps_point
{
  public:
    point a;    //点
    int index;  //点在原点集中的索引
};

记录每一个邻域点所属于的核心点的类：

class neps_list
{
  public:
    int c_idx;    //该邻域内点所属核心点的索引
    int n_idx;    //该邻域内点的索引
};

求曼哈顿距离代码：

float squareDistance(point a, point b)
{
  return ((abs(a.x-b.x)+abs(a.y-b.y))+abs(a.z-b.z));
}

将Mat矩阵载入point数组中的代码：

vector<point> load_mat_to_array(Mat img)
{
  vector<point> p;
  point tmp;
  for(int i = 0; i < img.rows; i++)
  {
    for(int j = 0; j < img.cols; j++)
    {
      tmp.x = (float)j;   //x坐标
      tmp.y = (float)i;   //y坐标
      tmp.z = (float)img.ptr<uchar>(i)[j];   //像素值
      p.push_back(tmp);
    }
  }
  return p;
}

把灰度图中的所有像素点分类之后，给每一类像素点进行着色的代码，其中每一类都着不同的颜色：

#define randomInt(a,b) (rand()%(b-a+1)+a)


void show_cluster_img(vector<vector<point>> Clusters, Mat &img, int row, int col)
{
  Mat img_tmp = Mat::zeros(row, col, CV_8UC3);


  srand((unsigned)time(NULL));
  for(int i = 0; i < Clusters.size(); i++)   //打印分类结果
  {
    int r = randomInt(0,255);
    int b = randomInt(0,255);
    int g = randomInt(0,255);
    for(int j = 0; j < Clusters[i].size(); j++)
    {
      int x = (int)Clusters[i][j].x;
      int y = (int)Clusters[i][j].y;
      
      if(x >= 0 && x < col && y >= 0 && y < row)
      {
        img_tmp.at<Vec3b>(y, x)[0] = r;
        img_tmp.at<Vec3b>(y, x)[1] = b;
        img_tmp.at<Vec3b>(y, x)[2] = g;
      }
    }
    
  }


  img_tmp.copyTo(img);


}

下面是主菜，DBSCAN算法代码：

vector<vector<point>> Dbscan(vector<point> &p, float Eps, int MinPts)
{
  vector<vector<neps_list>> c_p(p.size());   //记录每一个点的邻域点集
  for(int i=0; i < p.size(); i++)
  {
    for(int j=i; j < p.size(); j++)
    {
      if(squareDistance(p[i], p[j]) < Eps)
      {
        p[i].pts++;    //计数邻域内的点数
        neps_list t;
        t.c_idx = i;
        t.n_idx = j;
        c_p[i].push_back(t);   //将点j加入点i的邻域点集中


        if(i != j)   
        {
          p[j].pts++;
          t.c_idx = j;
          t.n_idx = i;
          c_p[j].push_back(t);    //将点i加入点j的邻域点集中
        }
      }
    }
  }


  
  for(int i = 0; i < p.size(); i++)   //判断邻域内的点数是否达到MinPts，达到则认为是核心点
  {
    if(p[i].pts >= MinPts)    //如果索引i的点是核心点，则标记其邻域所有点
    {
      p[i].pointType = 3;    //标记核心点
    }
  }


  vector<vector<point>> Clusters;    //簇集合


  int len = p.size();


  int cluster_num = -1;    //簇号初始化为-1


  for(int i = 0; i < len; i++)   //循环遍历每一个点
  {
    if(p[i].visited == 1)    //如果当前点已经被访问，则跳过
      continue;


    p[i].visited = 1;     //如果当前点未被访问，则标记为已访问


    if(p[i].pointType == 3)    //如果当前点为核心点
    {
      vector<point> C;   //新建一个簇
      cluster_num++;   //簇序号加1
      C.push_back(p[i]);   //将当前核心点加入到新建的簇中
      p[i].cluster = cluster_num;   //将当前的簇序号赋值给该点的所属簇序号
      
      vector<neps_point> N;  //求当前核心点的邻域点集合
      for(int k = 0; k < c_p[i].size(); k++)
      {
        neps_point tt;
        tt.a = p[c_p[i][k].n_idx];
        tt.index = c_p[i][k].n_idx;
        N.push_back(tt);
      }


      for(int j = 0; j < N.size(); j++)   //遍历邻域点集合中的所有点
      {
        if(p[N[j].index].visited == 0)   //通过index访问原点集中的当前邻域点，如果未被访问，则往下执行
        {
          p[N[j].index].visited = 1;   //在原点集中标记该点为已访问
          N[j].a.visited = 1;         //同时在邻域点集中标记该点为已访问


          if(p[N[j].index].pointType == 3)
          {
            for(int k = 0; k < c_p[N[j].index].size(); k++)
            {
              neps_point tt;
              tt.a = p[c_p[N[j].index][k].n_idx];
              tt.index = c_p[N[j].index][k].n_idx;
              N.push_back(tt);
            }
          }


          if(N[j].a.cluster == -1)   //如果当前遍历点尚未加入任何簇
          {
            C.push_back(p[N[j].index]);    //将该点加入新建的簇中
            N[j].a.cluster = cluster_num;   //将当前的簇序号赋值给该点的所属簇序号
            p[N[j].index].cluster = cluster_num;   //将当前的簇序号赋值给该点的所属簇序号
          }
        }
      }
      Clusters.push_back(C);    //将新建的簇加入簇集合中
      printf("Clusters.size()=%d\n", Clusters.size());
    }
 
  }


  return Clusters;
}

main函数代码：

int main()
{
  Mat img = imread("rgb.jpg", CV_LOAD_IMAGE_GRAYSCALE);
  resize(img, img, Size(round(img.cols*0.5), round(img.rows*0.5)), 2);
  imshow("img", img);


  vector<point> p = load_mat_to_array(img);


  printf("p.size()=%d\n", p.size());


  vector<vector<point>> Clusters = Dbscan(p, 8, 65);


  Mat rst;
  show_cluster_img(Clusters, rst, img.rows, img.cols);  //显示分类结果


  imshow("rst", rst);
  waitKey();


  return 0;
}

运行上述代码，得到的分类结果如下图所示。可以看到，图像被分成了8个区域，每个区域被涂上了不同的颜色，分类结果还是挺理想的。

实际使用该算法时，需要调整ε与minPts的值，以获取理想的分类结果，具体怎么调节，还没搞清楚，欢迎读者给我留言交流。