图像分割算法

阈值分割、连通域分析、区域生长、分水岭

阈值分割

固定阈值、自适应阈值(adaptiveThreshold)、大津阈值(OTSU)、最大熵阈值(KSW)

固定阈值

固定阈值的调用函数：

//InputArray类型的src，输入图像。
//OutputArray类型的dst，输出图像。
//double类型的thresh，阈值的具体值。
//double类型的maxval，阈值最大值。
//int类型的type，阈值操作的类型，包括：
//THRESH_BINARY（标准的二值化阈值法，大于thresh的设为maxval，小于的设为0），
//THRESH_BINARY_INV（反向二值化），
//THRESH_TRUNC（截断阈值法，大于thresh的设为thresh，小于则不变），
//THRESH_TOZERO（零化阈值法，大于thresh的不变，小于则零化），
//THRESH_TOZERO_INV（反向零化），
//THRESH_OTSU（大津算法，适合双峰直方图的图像，通过分析最大的背景前景类间方差，自动调节阈值），
//THRESH_TRIANGLE（三角法，适合单峰直方图图像，建立谷底和峰顶直线，距离直线垂直距离最大的直方图位置，即阈值thresh）。
double threshold( InputArray src, OutputArray dst, double thresh, double maxval, int type );

固定阈值的使用非常简单，例如：

cv::Mat im = cv::imread("1.jpg");
cv::Mat imBin;
cv::threshold(im, imBin, 80, 255, cv::THRESH_BINARY);

意思就是：im图像中大于80的像素值直接置为255，否则置为0，并保存在imBin中。

自适应阈值

简单来讲：就是用一个blockSize*blockSize大小的核，在图像中按像素单位移动，计算每次移动过程中的加权均值作为阈值的参考值。

//参数1:要处理的原图
//参数2: 最大闯值，一般为255
//参数3: 小区域闽值的计算方式:
//ADAPTIVE THRESH MEAN C: 小区域内取均值
//ADAPTIVE THRESH GAUSSIAN C: 小区域内加权求和，权重是个高斯核
//参数4: 这是阈值类型，只有两个取值，分别为 THRESH_BINARY 和THRESH_BINARY_INV。
//参数5: 小区域的面积，如11就是11*11的小块
//参数6: 最终闻值等于小区域计算出的闯值再减去此值
//特定: 自适应闽值会每次取图片的一小部分计算闽值，这样图片不同区域的闽值就不尽相同，适用于明暗分布不均的图片。
void adaptiveThreshold(InputArray src, OutputArray dst, double maxValue, int adaptiveMethod, int thresholdType, int blockSize, double C)

主要是用于光照变化明显的图像，但针对于一般场景的图像分割，效果往往不甚理想。
简单使用：

cv::Mat im = cv::imread("1.jpg");
cv::Mat imBin;
cv::adaptiveThreshold(im, imBin, 255, cv::ADAPTIVE_THRESH_MEAN_C, cv::THRESH_BINARY, 21, 10);

稍加修改：

enum adaptiveMethod{meanFilter,gaaussianFilter,medianFilter};
void AdaptiveThreshold(cv::Mat& src, cv::Mat& dst, double Maxval, int Subsize, double c, adaptiveMethod method = meanFilter){if (src.channels() > 1)cv::cvtColor(src, src, CV_RGB2GRAY);cv::Mat smooth;switch (method){case  meanFilter:cv::blur(src, smooth, cv::Size(Subsize, Subsize));  //均值滤波break;case gaaussianFilter:cv::GaussianBlur(src, smooth, cv::Size(Subsize, Subsize),0,0); //高斯滤波break;case medianFilter:cv::medianBlur(src, smooth, Subsize);   //中值滤波break;default:break;} smooth = smooth - c;//阈值处理src.copyTo(dst);for (int r = 0; r < src.rows;++r){const uchar* srcptr = src.ptr<uchar>(r);const uchar* smoothptr = smooth.ptr<uchar>(r);uchar* dstptr = dst.ptr<uchar>(r);for (int c = 0; c < src.cols; ++c){if (srcptr[c]>smoothptr[c]){dstptr[c] = Maxval;}elsedstptr[c] = 0;}}
}

大津阈值(OTSU)

简单来讲就是利用类间方差的方式，找到背景和目标差异最大的阈值。因方差是灰度分布均匀性的一种度量,背景和前景之间的类间方差越大,说明构成图像的两部分的差别越大,当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。
OpenCV中也包含了大津阈值的的使用接口：

cv::Mat im = cv::imread("1.jpg");
cv::Mat imBin;
double thresh=0;
thresh = cv::threshold(im, imBin, thresh, 255, CV_THRESH_OTSU + CV_THRESH_BINARY);

针对于一般的图像分割时，不会有太大的问题，但是如果是背景太多的图像，如下图所示：

采用上述方式得到的结果为：

不能将图像中的裂缝分割出来。其实那些黑色背景并不需要添加到大津阈值的计算当中，因此我们需要了解大津阈值的具体计算过程，并修改源代码。
OTSU算法的假设是存在阈值K将图像所有像素分为前景$N_1$(小于K)背景$N_2$(大于K)。假设图像大小为M*N，大于阈值K的像素为$I_1(x,y)$，小于阈值K的像素为$I_2(x,y)$则这两类像素各自的均值就为$m_1$、$m_2$：
$$m_1=\frac{\sum I_1(x,y)}{N_1};m_2=\frac{\sum I_2(x,y)}{N_2}$$属于前景的概率$P_1$和属于背景的概率$P_2$分别为：
$$P_1=\frac{N_1}{\mathrm{M}\ast \mathrm{N}};P_2=\frac{N_2}{\mathrm{M}\ast \mathrm{N}}$$且：
$$P_1+P_1=1$$图像全局均值为$m_G$：
$$m_G=\frac{\sum I(x,y)}{M\ast N}$$且：
$$P_1{m}_1+P_1{m}_2=m_G$$根据方差的概念，类间方差表达式为：
$${\sigma }^2=P_1(m_1-m_G{)}^2+P_2(m_2-m_G{)}^2$$背景和目标之间的类间方差越大，说明构成图像的两部分的差别越大，当部分目标错分为背景或部分背景错分为目标都会导致两部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。
将上式稍微变形可以得到：
$${\sigma }^2=P_1{P}_2(m_1-m_2{)}^2$$具体推到过程如下图所示：

照着公式，遍历0~255个灰度级，求出使上式类间方差最大灰度级作为 K值就可以了。但还可以对上式进行如下变化，使得式子更加简单，推到过程如下：
假设$p_i$表示每个灰度级的概率(i=0~255)，$N_i$表示每个灰度级的数量，$N_{1i}$表示前景部分每个灰度级的数量，$N_{2i}$表示背景部分每个灰度级的数量。
$$p_i=\frac{N_i}{M\ast N}$$因此可以得到：
$$P_1=\sum _{i=0}^K{p}_i,P_2=\sum _{i=K+1}^{255}{p}_i$$$$m_1=1/P_1\ast \sum _{i=0}^{K}i{p}_i=\frac{\sum I_1(x,y)}{N_1}$$$$m_2=1/P_2\ast \sum _{i=K+1}^{255}i{p}_i=\frac{\sum I_2(x,y)}{N_2}$$定义灰度级为K时前景的累加均值m为：
$$m=\sum _{i=0}^{K}i{p}_i$$图像全局均值$m_G$又可以写为：$$m_G=\sum _{i=0}^{255}i{p}_i=\frac{\sum I(x,y)}{M\ast N}$$因此可以改写$m_1$和$m_2$为：
$$m_1=1/P_1\ast m$$$$m_2=1/P_2\ast (m_G-m)$$r然后可以将类间方差公式改写为：
$${\sigma }^2=\frac{{(m_G\ast P_1-m)}^2}{P_1(1-P_1)}$$具体推到过程如下图所示：

最后贴上代码

int Otsu(cv::Mat& src, cv::Mat& dst)
{const int Grayscale = 256;int graynum[Grayscale] = { 0 };int r = src.rows;int c = src.cols;double srcpixnum = 0;for (int i = 0; i < r; ++i){uchar* ptr = src.ptr<uchar>(i);for (int j = 0; j < c; ++j) {        //直方图统计if (ptr[j] == 0)  //不计算背景{ptr[j] = uchar(255); //背景变亮(根据实际情况而定)continue;}graynum[ptr[j]]++;srcpixnum++;}}double P[Grayscale] = { 0 };double PK[Grayscale] = { 0 };double MK[Grayscale] = { 0 };double sumtmpPK = 0, sumtmpMK = 0;for (int i = 0; i < Grayscale; ++i) {P[i] = graynum[i] / srcpixnum;   //每个灰度级出现的概率PK[i] = sumtmpPK + P[i];         //概率累计和 sumtmpPK = PK[i];MK[i] = sumtmpMK + i * P[i];       //灰度级的累加均值                                                                                                                                                                                                                                                                                                                                                                                                        sumtmpMK = MK[i];}//计算类间方差int thresh;double Var = 0;for (int k = 0; k < Grayscale; ++k) {if ((MK[Grayscale - 1] * PK[k] - MK[k])*(MK[Grayscale - 1] * PK[k] - MK[k]) / (PK[k] * (1 - PK[k])) > Var) {Var = (MK[Grayscale - 1] * PK[k] - MK[k])*(MK[Grayscale - 1] * PK[k] - MK[k]) / (PK[k] * (1 - PK[k]));thresh = k;}}//阈值处理dst = src.clone();for (int i = 0; i < r; ++i) {uchar* ptr = dst.ptr<uchar>(i);for (int j = 0; j < c; ++j){if (ptr[j] > thresh)  //大于阈值变暗(根据实际情况而定)ptr[j] = 0;elseptr[j] = 255;}}return thresh;
}

根据上述代码进行上面图像分割可以得到：

最大熵阈值

最大熵阈值分割法和OTSU算法类似，假设将图像分为背景和前景两个部分。熵代表信息量，图像信息量越大，熵就越大，最大熵算法就是找出一个最佳阈值使得背景与前景两个部分熵之和最大。

连通域分析

连通区域(Connected Component) 一般是指图像中具有相同像素值且位置相邻的前景像素点组成的图像区域，连通区域分析是指将图像中的各个连通区域找出并标记。在图像中，最小的单位是像素，每个像素周围有邻接像素，常见的连通关系有2种: 4连通与8连通。


如果A与B连通，B与C连通，则A与C连通，在视觉上看来，彼此连通的点形成了一个区域，而不连通的点形成了不同的区域。这样的一个所有的点彼此连通点构成的集合，我们称为一个连通区域。

4邻接，则有3个连通区域；8邻接，则有2个连通区域。

OpenCV的连通与分析，参考上面链接。这里简单介绍一下Two-Pass算法。两遍扫描法( Two-Pass )，正如其名，指的就是通过扫描两遍图像，将图像中存在的所有连通域找出并标记。

第一次扫描 :
从左上角开始遍历像素点，找到第一个像素为255的点label=1 ；
当该像素的左邻像素和上邻像素为无效值时，给该像素置个新的label值，label ++，记录集合；
当该像素的左邻像素或者上邻像素有一个为有效值时，将有效值像素的label赋给该像素的label值 ；
当该像素的左邻像素和上邻像素都为有效值时，选取其中较小的label值赋给该像素的label值。

第二次扫描 :
对每个点的label进行更新，更新为其对于其集合中最小的label。

区域生长

基本思想是将具有相似性质的像素集合起来构成区域步骤 :

1. 在图像中选取若干个种子点，将种子点存入vector中，vector中存储待生长的种子点；选取一个种子点(x0,y0);
2. 以(x0,y0)为中心,考虑(x0,y0)的4邻域像素(x, y)如果(x0,y0)满足生长准则,将(x, y)与(x0,y0)合并(在同一区域内)，同时将(x,y)压入堆栈;
3. 从堆栈中取出一个像素,把它当作(x0,y0)返回到步骤2:
4. 当堆栈为空时，返回到步骤1;
5. 重复步骤1 - 4直到图像中的每个点都有归属时;
6. 生长结束。

/***************************************************************************************
Function:  区域生长算法
Input:     src 待处理原图像 pt 初始生长点 th 生长的阈值条件
Output:    肺实质的所在的区域 实质区是白色，其他区域是黑色
Description: 生长结果区域标记为白色(255),背景色为黑色(0)
Return:    NULL
Others:    NULL
***************************************************************************************/
void RegionGrow(cv::Mat& src, cv::Mat& matDst, cv::Point2i pt, int th)
{cv::Point2i ptGrowing;							//待生长点位置int nGrowLable = 0;								//标记是否生长过int nSrcValue = 0;								//生长起点灰度值int nCurValue = 0;								//当前生长点灰度值matDst = cv::Mat::zeros(src.size(), CV_8UC1);	//创建一个空白区域，填充为黑色//生长方向顺序数据int DIR[8][2] = { { -1, -1 }, { 0, -1 }, { 1, -1 }, { 1, 0 }, { 1, 1 }, { 0, 1 }, { -1, 1 }, { -1, 0 } };std::vector<cv::Point2i> vcGrowPt;				//生长点栈vcGrowPt.push_back(pt);							//将生长点压入栈中matDst.at<uchar>(pt.y, pt.x) = 255;				//标记生长点nSrcValue = src.at<uchar>(pt.y, pt.x);			//记录生长点的灰度值while (!vcGrowPt.empty())						//生长栈不为空则生长{pt = vcGrowPt.back();						//取出一个生长点vcGrowPt.pop_back();//分别对八个方向上的点进行生长for (int i = 0; i < 8; ++i){ptGrowing.x = pt.x + DIR[i][0];ptGrowing.y = pt.y + DIR[i][1];//检查是否是边缘点if (ptGrowing.x < 0 || ptGrowing.y < 0 || ptGrowing.x >(src.cols - 1) || (ptGrowing.y > src.rows - 1))continue;nGrowLable = matDst.at<uchar>(ptGrowing.y, ptGrowing.x);		//当前待生长点的灰度值if (nGrowLable == 0)					//如果标记点还没有被生长{nCurValue = src.at<uchar>(ptGrowing.y, ptGrowing.x);if (abs(nSrcValue - nCurValue) < th)					//在阈值范围内则生长{matDst.at<uchar>(ptGrowing.y, ptGrowing.x) = 255;		//标记为白色vcGrowPt.push_back(ptGrowing);					//将下一个生长点压入栈中}}}}
}

分水岭

分水岭算法步骤：
1、加载原始图像
2、闽值分割，将图像分割为黑白两个部分
3、对图像进行开运算，即先腐蚀在膨胀
4、对开运算的结果再进行 膨胀，得到大部分是背景的区域
5、通过距离变换 Distance Transform 获取前景区域
6、背景区域sure bg 和前景区域sure fg相减，得到即有前景又有背景的重合区域
7、连通区域处理
8、最后使用分水岭算法
一个使用分水岭算法的例子：

//分水岭算法
Mat WaterSegment(Mat src) {int row = src.rows;int col = src.cols;//1. 将RGB图像灰度化Mat grayImage;cv::cvtColor(src, grayImage, cv::COLOR_BGR2GRAY);//2. 使用大津法转为二值图，并做形态学闭合操作threshold(grayImage, grayImage, 0, 255, THRESH_BINARY | THRESH_OTSU);//3. 形态学闭操作Mat kernel = getStructuringElement(MORPH_RECT, Size(9, 9), Point(-1, -1));morphologyEx(grayImage, grayImage, MORPH_CLOSE, kernel);//4. 距离变换distanceTransform(grayImage, grayImage, DIST_L2, DIST_MASK_3, 5);//5. 将图像归一化到[0, 1]范围normalize(grayImage, grayImage, 0, 1, NORM_MINMAX);//6. 将图像取值范围变为8位(0-255)grayImage.convertTo(grayImage, CV_8UC1);//7. 再使用大津法转为二值图，并做形态学闭合操作threshold(grayImage, grayImage, 0, 255, THRESH_BINARY | THRESH_OTSU);morphologyEx(grayImage, grayImage, MORPH_CLOSE, kernel);//8. 使用findContours寻找marksvector<vector<Point>> contours;findContours(grayImage, contours, RETR_TREE, CHAIN_APPROX_SIMPLE, Point(-1, -1));Mat marks = Mat::zeros(grayImage.size(), CV_32SC1);for (size_t i = 0; i < contours.size(); i++){//static_cast<int>(i+1)是为了分水岭的标记不同，区域1、2、3...这样才能分割drawContours(marks, contours, static_cast<int>(i), Scalar::all(static_cast<int>(i + 1)), 2);}//9. 对原图做形态学的腐蚀操作Mat k = getStructuringElement(MORPH_RECT, Size(3, 3), Point(-1, -1));morphologyEx(src, src, MORPH_ERODE, k);//10. 调用opencv的分水岭算法watershed(src, marks);//11. 随机分配颜色vector<Vec3b> colors;for (size_t i = 0; i < contours.size(); i++) {int r = theRNG().uniform(0, 255);int g = theRNG().uniform(0, 255);int b = theRNG().uniform(0, 255);colors.push_back(Vec3b((uchar)b, (uchar)g, (uchar)r));}// 12. 显示Mat dst = Mat::zeros(marks.size(), CV_8UC3);int index = 0;for (int i = 0; i < row; i++) {for (int j = 0; j < col; j++) {index = marks.at<int>(i, j);if (index > 0 && index <= contours.size()) {dst.at<Vec3b>(i, j) = colors[index - 1];}else if (index == -1){dst.at<Vec3b>(i, j) = Vec3b(255, 255, 255);}else {dst.at<Vec3b>(i, j) = Vec3b(0, 0, 0);}}}return dst;
}

总结：个人觉得分水岭算法其实很依赖图像阈值分割算法的处理效果，然后就是轮廓检测，这两步没处理好的话，后续效果也一般。
参考：https://blog.csdn.net/just_sort/article/details/87376355