【图像处理笔记】图像分割之聚类和超像素

0 引言

大多数分割算法都基于图像灰度值的两个基本性质之一:不连续性相似性第一类方法根据灰度的突变(如边缘)将图像分割为多个区域:首先寻找边缘线段,然后将这些线段连接为边界的方法来识别区域。第二类方法根据一组预定义的准则把一幅图像分割为多个区域。本节讨论两种相关的区域分割方法:(1)在数据中寻找聚类的经典方法,它与亮度和颜色等变量有关;(2)用聚类从图像中提取“超像素”的现代方法。

【图像处理笔记】图像分割之聚类和超像素插图

1 使用k均值聚类的区域分割

1.1 原理

聚类方法的思想是将样本集合按照其特征的相似性划分为若干类别,使同一类别样本的特征具有较高的相似性,不同类别样本的特征具有较大的差异性。令{z1, z2, z3 ..., zn}是样本集合,在图像分割中,样本向量z的每个分量表示一个数值像素属性。例如,分割只基于灰度尺度时,z是一个表示像素灰度的标量。分割的如果是RGB彩色图像,z通常是一个三维向量,这个三维向量的每个分量是RGB三通道的灰度值。k均值聚类的目的就是将样本集合划分为k个满足如下最优准则的不相交的聚类集合C={C1, C2, ..., Ck}:

【图像处理笔记】图像分割之聚类和超像素插图1

式中,mi是集合Ci中样本的均值向量(或质心),||z-mi||项是Ci中的一个样本到均值mi的欧式距离。换言之,这个公式说,我们感兴趣的是找到集合C={C1, C2, ..., Ck},集合中的每个点到该集合的均值的距离之和是最小的。

基于聚类的区域分割,就是基于图像的灰度颜色纹理形状等特征,用聚类算法把图像分成若干类别或区域,使每个点到聚类中心的均值最小。k 均值(k-means)是一种无监督聚类算法。基于 k 均值聚类算法的区域分割,算法步骤为:

(1)初始化算法:规定一组初始均值

(2)将样本分配给聚类:对所有的像素点,计算像素到每个聚类中心的距离,将像素分类到距离最小的一个聚类中;

(3)更新聚类中心:根据分类结果计算出新的聚类中心;

(4)完备性验证:计算当前步骤和前几步中平均向量之间的差的欧几里得范数。计算残差E,即k个范数之和。若E≤T,其中T是一个规定的非负阈值,则停止。否则,返回步骤2。

1.2 cv::kmeans函数

OpenCV提供了函数cv::kmeans来实现 k-means 聚类算法。函数cv::kmeans不仅可以基于灰度、颜色对图像进行区域分割,也可以基于样本的其它特征如纹理、形状进行聚类。

double cv::kmeans(InputArray data,  //用于聚类的数据,类型为 CV_32F
                  int K, //设定的聚类数量
                  InputOutputArray bestLabels,  //输出整数数组,用于存储每个样本的聚类类别索引
                  TermCriteria criteria,  //算法终止条件:即最大迭代次数或所需精度
                  int attempts,  //用于指定使用不同初始标记执行算法的次数
                  int flags,  //初始化聚类中心的方法:0=随机初始化 1=kmeans++方法初始化 2=第一次用用户指定的标签初始化,后面attempts-1都用随机或版随机的初始化
                  OutputArray centers = noArray()  //聚类中心的输出矩阵,每个聚类中心占一行
                  )

示例 图像分割之k均值聚类

Mat src = imread("./14.tif", 0);
Mat dataPixels = src.reshape(0, 1);//可以是一列,每一行表示一个样本;或者一行,一列是一个样本;样本的分量数为通道数
dataPixels.convertTo(dataPixels, CV_32FC1);//输入需要是32位浮点型
int numCluster = 3;
Mat labels;
Mat centers;
kmeans(dataPixels, numCluster, labels, TermCriteria(TermCriteria::EPS + TermCriteria::COUNT, 10, 0.1), 3, KMEANS_PP_CENTERS, centers);
Mat dst = Mat::zeros(src.size(), CV_8UC1);
float* pdata = dataPixels.ptr(0);
for (int i = 0; i (i)[0];//每个像素对应的标签k,即属于集合k
    pdata[i] = centers.ptr(k)[0];//用集合中心替换该像素
}
dataPixels.convertTo(dataPixels, CV_8UC1);
dst = dataPixels.reshape(0, src.rows);

【图像处理笔记】图像分割之聚类和超像素插图2

1.3 cv::kmeans源码

当图片非常大时,对图像进行简单的计算操作,耗时就会变得非常大,常用的加速方法如OpenMp,TBB,OpenCL和CUDA。使用cuda时,图片在cpu和gpu之间的传输时间就达到上百ms,不适合本来就是ms级的计算。在kmeans源码中使用parallel_for_并行计算各个样本到聚类中心的距离。这边写一个简单的例子,了解下parallel_for_的用法。

示例 利用并行计算加速图片旋转

// 方法1:并行 将该方法写成一个类,继承ParallelLoopBody,然后重写(),利用parallel_for_可以开启并行
class trans :public ParallelLoopBody {
public:
    trans(const uchar* _src, uchar*_dst, int _dims, int _istep):src(_src), dst(_dst), dims(_dims), istep(_istep) {}
    void operator()(const Range& range) const //重载操作符()
    {
        for (int n = range.start; n ();
    uchar* pdst = dst.ptr();
    for (size_t i = 0; i (); 
uchar* pdst = dst.ptr();
int dims = src.cols;
int istep = src.step;
int N = src.rows;
rotate(src, dst, ROTATE_90_CLOCKWISE);// opencv提供的方法 3ms
rotate(src, dst, src.cols, src.rows);// 自己写的for循环 20ms
parallel_for_(Range(0, N), trans(psrc, pdst, dims, istep)); // 并行加速3ms

源码位于opencv路径下sourcesmodulescoresrckmeans.cpp中,首先是初始化算法:规定一组初始center,一种是随机产生,另一种是用kmeans++初始化。kmeans++初始化聚类中心是以概率的形式逐个选择聚类中心,并在选择聚类中心时,给距离较远的点更高的权重,即更容易被选择为聚类中心。假设有5个点,随机选择其中一个点为中心,计算其他点到该点的距离的平方分别为10,20,5,15。则选下一个聚类中心时它们的权值为0.2,0.4,0.1和0.3。用代码写就是,距离平方和50,在【0,50】间随机生成一个数,用这个数挨个减去10,再减去20...直到结果小于0,下一个聚类中心就是该点。

// 随机产生聚类中心
// dims 样本向量的分量数
// box  存放了所有样本中每个分量的最大值和最小值
static void generateRandomCenter(int dims, const Vec2f* box, float* center, RNG& rng)
{
    float margin = 1.f / dims;
    for (int j = 0; j  _centers(K);
    int* centers = &_centers[0];
    cv::AutoBuffer _dist(N * 3);// 3倍样本数量大小,存不同时刻样本到最近聚类中心的距离
    float* dist = &_dist[0], * tdist = dist + N, * tdist2 = tdist + N;
    double sum0 = 0;
    // 第一个聚类中心随机生成
    centers[0] = (unsigned)rng % N;// %N 即可获得[0,N-1]的随机数
    // 计算所有样本到第一个中心的距离,并求和
    for (int i = 0; i (i), data.ptr(centers[0]), dims);
        sum0 += dist[i];
    }
    // 计算第2、3...个聚类中心,离第一个中心越远的点权重越高
    for (int k = 1; k (centers[k]);
        float* dst = _out_centers.ptr(k);
        for (int j = 0; j 

计算距离

// 并行计算每个样本距离中心的距离
// dist 3倍样本数量N的autobuffer,前N个存放上次的N个距离(dist_指向第一个),后N个存放本次计算的N个距离(tdist2_指向第一个)
// data 样本向量集合
// ci   计算所有样本和第ci个样本的距离
class KMeansPPDistanceComputer : public ParallelLoopBody
{
public:
    KMeansPPDistanceComputer(float* tdist2_, const Mat& data_, const float* dist_, int ci_) :
        tdist2(tdist2_), data(data_), dist(dist_), ci(ci_)//成员初始化列表的写法
    { }

    void operator()(const cv::Range& range) const CV_OVERRIDE
    {
        CV_TRACE_FUNCTION();
        const int begin = range.start;
        const int end = range.end;
        const int dims = data.cols;

        for (int i = begin; i (i), data.ptr(ci), dims), dist[i]);
        }
    }

private:
    KMeansPPDistanceComputer& operator=(const KMeansPPDistanceComputer&); // = delete

    float* tdist2;
    const Mat& data;
    const float* dist;
    const int ci;
};

更新标签

// 并行计算每个样本到对应中心的距离,已知样本属于哪个集合,直接计算该样本到集合中心的距离
template
class KMeansDistanceComputer : public ParallelLoopBody
{
public:
    KMeansDistanceComputer(double* distances_,
        int* labels_,
        const Mat& data_,
        const Mat& centers_)
        : distances(distances_),
        labels(labels_),
        data(data_),
        centers(centers_)
    {
    }

    void operator()(const Range& range) const CV_OVERRIDE
    {
        CV_TRACE_FUNCTION();
        const int begin = range.start;
        const int end = range.end;
        const int K = centers.rows;
        const int dims = centers.cols;

        for (int i = begin; i (i);
            if (onlyDistance)// 只算距离,不重新分配标签
            {
                const float* center = centers.ptr(labels[i]);
                distances[i] = hal::normL2Sqr_(sample, center, dims);
                continue;
            }
            else // 遍历每一个样本,计算该样本到每一个中心的距离,重新分配标签
            {
                int k_best = 0;
                double min_dist = DBL_MAX;

                for (int k = 0; k (k);
                    const double dist = hal::normL2Sqr_(sample, center, dims);

                    if (min_dist > dist)
                    {
                        min_dist = dist;
                        k_best = k;
                    }
                }

                distances[i] = min_dist;
                labels[i] = k_best;
            }
        }
    }

private:
    KMeansDistanceComputer& operator=(const KMeansDistanceComputer&); // = delete

    double* distances;
    int* labels;
    const Mat& data;
    const Mat& centers;
};

kmeans

//_data:特征向量集;K:聚类中心个数;_bestLabels:每个特征向量隶属聚类中心索引
 //criteria:迭代算法终止条件;attempts:尝试次数;
 //flags:第一次尝试初始化采取策略;_centers:存放聚类中心
double cv::kmeans(InputArray _data, int K,
    InputOutputArray _bestLabels,
    TermCriteria criteria, int attempts,
    int flags, OutputArray _centers)
{
    CV_INSTRUMENT_REGION();
    const int SPP_TRIALS = 3;
    Mat data0 = _data.getMat();
    const bool isrow = data0.rows == 1;// 输入的数据应该是一行,或者一列的,通道数是每个样本向量的分量数
    const int N = isrow ? data0.cols : data0.rows;// N表示样本向量个数
    const int dims = (isrow ? 1 : data0.cols) * data0.channels();// 每个样本向量的维度(分量数)
    const int type = data0.depth();//数据类型,应为32位浮点数

    attempts = std::max(attempts, 1);//至少尝试一次
    CV_Assert(data0.dims  0);
    CV_CheckGE(N, K, "Number of clusters should be more than number of elements");

    Mat data(N, dims, CV_32F, data0.ptr(), isrow ? dims * sizeof(float) : static_cast(data0.step));

    _bestLabels.create(N, 1, CV_32S, -1, true);//_bestLabels为N*1大小矩阵,类型为32为有符号整型,每个样本向量有有一个标签

    Mat _labels, best_labels = _bestLabels.getMat();
    if (flags & CV_KMEANS_USE_INITIAL_LABELS) // 如果首次是由用户指定的
    {
        CV_Assert((best_labels.cols == 1 || best_labels.rows == 1) &&
            best_labels.cols * best_labels.rows == N &&
            best_labels.type() == CV_32S &&
            best_labels.isContinuous());
        best_labels.reshape(1, N).copyTo(_labels);
        for (int i = 0; i (i) ();

    Mat centers(K, dims, type), old_centers(K, dims, type), temp(1, dims, type);
    cv::AutoBuffer counters(K);
    cv::AutoBuffer dists(N);
    RNG& rng = theRNG();

    if (criteria.type & TermCriteria::EPS)//算法精度
        criteria.epsilon = std::max(criteria.epsilon, 0.);
    else
        criteria.epsilon = FLT_EPSILON;
    criteria.epsilon *= criteria.epsilon;

    if (criteria.type & TermCriteria::COUNT)//最大迭代次数
        criteria.maxCount = std::min(std::max(criteria.maxCount, 2), 100);
    else
        criteria.maxCount = 100;

    if (K == 1)
    {
        attempts = 1;
        criteria.maxCount = 2;
    }

    cv::AutoBuffer box(dims);
    if (!(flags & KMEANS_PP_CENTERS))//随机初始化中心的话要计算下范围,在最大小值之间随机生成
    {
        {
            const float* sample = data.ptr(0);
            for (int j = 0; j (i);
            for (int j = 0; j  0 || !(flags & KMEANS_USE_INITIAL_LABELS)))
            {
                if (flags & KMEANS_PP_CENTERS)
                    generateCentersPP(data, centers, K, rng, SPP_TRIALS);
                else
                {
                    for (int k = 0; k (k), rng);
                }
            }
            else //若为人工指定labels,或者不是第一次迭代,将样本划分进不同的集合,根据labels
            {
                // compute centers
                centers = Scalar(0); // 对centers进行初始化操作
                for (int k = 0; k (i);
                    int k = labels[i];
                    float* center = centers.ptr(k);
                    for (int j = 0; j (max_k);
                    float* _base_center = temp.ptr(); // normalized
                    float scale = 1.f / counters[max_k];
                    for (int j = 0; j (i);
                        double dist = hal::normL2Sqr_(sample, _base_center, dims);

                        if (max_dist (farthest_i);
                    float* cur_center = centers.ptr(k);
                    for (int j = 0; j (k);
                    CV_Assert(counters[k] != 0);

                    float scale = 1.f / counters[k];
                    for (int j = 0; j  0)// 计算本次循环和上次聚类中心的差距,差距小于criteria.epsilon则为最后一次迭代
                    {
                        double dist = 0;
                        const float* old_center = old_centers.ptr(k);
                        for (int j = 0; j (dists.data(), labels, data, centers), (double)divUp((size_t)(dims * N), CV_KMEANS_PARALLEL_GRANULARITY));
                compactness = sum(Mat(Size(N, 1), CV_64F, &dists[0]))[0];// 记录距离和
                break;
            }
            else // 不是最后一次的话,计算距离的同时还要重新分配下标签,可能会导致空集合
            {
                // assign labels
                parallel_for_(Range(0, N), KMeansDistanceComputer(dists.data(), labels, data, centers), (double)divUp((size_t)(dims * N * K), CV_KMEANS_PARALLEL_GRANULARITY));
            }
        }
        //compactness将记录所有距离,这里的距离是指,所有的特征向量到其聚类中心的距离之和,用于评价当前的聚类结果
        if (compactness 

 

2 使用超像素的区域分割

超像素图像分割基于依赖于图像的颜色信息空间关系信息将图像分割为远超于目标个数、远小于像素数量的超像素块,达到尽可能保留图像中所有目标的边缘信息的目的,从而更好的辅助后续视觉任务(如目标检测、目标跟踪、语义分割等)。

超像素是由一系列位置相邻,颜色、亮度、纹理等特征相似的像素点组成的小区域,我们将其视为具有代表性的大“像素”,称为超像素。超像素技术通过像素的组合得到少量(相对于像素数量)具有感知意义的超像素区域,代替大量原始像素表达图像特征,可以极大地降低图像处理的复杂度、减小计算量。超像素分割的结果是覆盖整个图像的子区域的集合,或从图像中提取的轮廓线的集合。 超像素的数量越少,丧失的细节特征越多,但仍然能基本保留主要区域之间的边界及图像的基本拓扑

【图像处理笔记】图像分割之聚类和超像素插图3

常用的超像素分割方法有:简单线性迭代聚类(Simple Linear Iterative Clustering,SLIC)、能量驱动采样(Super-pixels Extracted via Energy-Driven Sampling,SEEDS)和线性谱聚类(Linear Spectral Clustering,LSC)。SLIC超像素算法是对上节讨论的k均值算法的一种改进,通常使用(但不限于)包含三个颜色分量和两个空间坐标的五维向量。

OpenCV 在 ximgproc 模块提供了ximgproc.createSuperpixelSLIC函数实现SLIC算法。需要编译opencv_contrib模块,可以参考VS2019编译Opencv4.6.0GPU版本,记得勾选ximgproc。

示例 SLIC超像素区域分割

#include
using namespace ximgproc;
...
Mat src = imread("./14.tif");
Mat slicLabel, slicMask, slicColor, slicDst;
Ptr slic = createSuperpixelLSC(src);
slic->iterate(10);//迭代次数
slic->getLabels(slicLabel);//获取labels
slic->getLabelContourMask(slicMask);//获取超像素的边界
int number = slic->getNumberOfSuperpixels();//获取超像素的数量
src.setTo(Scalar(255, 255, 255), slicMask);

【图像处理笔记】图像分割之聚类和超像素插图4

 

 

参考 

1. 冈萨雷斯《数字图像处理(第四版)》Chapter 10(所有图片可在链接中下载)

2. 【youcans 的 OpenCV 例程200篇】171.SLIC 超像素区域分割

3. 当我们在谈论kmeans(3)

 

文章来源于互联网:【图像处理笔记】图像分割之聚类和超像素

THE END
分享
二维码