【图像处理笔记】图像分割之阈值处理
本章的大多数分割算法都基于图像灰度值的两个基本性质之一：不连续性和相似性。第一类方法根据灰度的突变（如边缘）将图像分割为多个区域；第二类方法根据一组预定义的准则把一幅图像分割为多个区域。上一节采用第一类方法首先寻找边缘线段，然后将这些线段连接为边界的方法来识别区域。本节讨论根据灰度值或灰度值的性质来将图像直接划分为多个区域的技术。由于图像阈值处理直观、实现简单并且计算速度快，因此在图像分割应用中处于核心地位。
1 基础知识
1.1 灰度阈值处理基础
假设下图(a)中的灰度直方图对应于图像f(x,y)，将图像由暗色背景上的亮目标组成，其中背景目标和目标像素组成了两种主导模式。从背景中提取目标的一种明显方法是，选择一个分隔这些模式的阈值T。然后，图像中f(x,y)大于T的任何点(x,y)称为一个目标点；否则该点称为背景点。换句话说，分割后的图像g(x,y)为
当T是一个适合于整个图像的常数时，上式中给出的处理称为全局阈值处理。当T值在一幅图像上变化时，上式中给出的处理称为可变阈值处理。有时，我们使用局部阈值处理或区域阈值处理表示可变阈值处理，此时，图像中任意一点(x,y)处的T值取决于(x,y)的邻域的性质（如邻域中像素的平均灰度）。若T取决于空间坐标(x,y)本身，则可变阈值处理通常称为动态阈值处理或自适应阈值处理。
上图(b)显示了一个更难的阈值处理问题——双阈值处理，它包含一个具有3个主导模式的直方图，这3个主导模式分别对应于暗色背景和两类亮目标。分割后的图像为
根据前面的讨论，我们可以凭直觉推断灰度阈值处理成功与否，与分隔直方图模式的波谷的宽度和深度有关。而影响波谷特性的关键因素依次是：（1）波峰间的间隔（波峰离得越远，分隔这些模式的机会越好）；（2）图像中的噪声内容（噪声越大，模式越宽）；（3）目标和背景的相对大小；（4）光源的均匀性；（5）图像反射性质的均匀性。
1.2 噪声的影响
下图(a)中没有噪声，因此其直方图由两个波峰模式组成，如图(d)所示。将该图像分割成两个区域很简单：我们只需要在两个模式之间的任何位置选择一个阈值。图(b)显示了被均值为零、标准差为10个灰度级的高斯噪声污染的原图像。现在，模式较宽(如图(e)所示)，但它们的间隔足以让它们之间的波谷的深度轻易的分隔两个模式。放在两个波峰之间的中间位置的一个阈值就能很好地分割该图像。图(c)显示了该图像被均值为零、标准差为50个灰度级的高斯噪声污染的结果。如图(f)所示，现在的情况严重到无法区分两个模式，找不到合适的阈值分割图像。
1.3 光照和反射的影响
下图(a)是一个带噪声图像，使用单个阈值很容易分割这幅图像。假设(a)中的图像乘以一个不均匀的灰度函数，如图(b)中的灰度斜坡，其直方图如(e)所示。图(c)显示了这两幅图的乘积，图(f)是得到的直方图。如果不进行额外的处理，那么波峰之间的深谷会破坏到模式无法分隔的程度。如果光照非常均匀，但（目标和/或背景表面自然反射率的变换导致的）图像反射不均匀，那么也会得到类似的结果。
照明和反射引起的不均匀性，可以用如下三种基本方法解决：
1. 直接校正这种阴影模式。例如，不均匀（但固定的）光照可以用相反的模式与图像相乘来校正，相反的模式可以通过对以一个恒定灰度的平坦表面成像来得到。
2. 用顶帽变换处理来校正全局阴影模式。
3. 用可变阈值处理来“绕过”不均匀性
2 基本的全局阈值处理
当目标和背景像素的灰度分布非常不同时，可对整个图像使用单个（全局）阈值。在大多数应用中，图像之间通常存在足够的变化，即使全局阈值是一种合适的方法，也需要有能对每幅图像估计阈值的算法。下面的迭代算法适用于这一目的：
（1）设置初始阈值 T，通常可以设为图像的平均灰度；
（2）用灰度阈值 T 分割图像：灰度值等于 T 的所有像素集合 G1 和 大于等于 T 的所有像素集合 G2；
（3）分别计算 G1、G2 的平均灰度值 m1、m2；
（4）求出新的灰度阈值 T=(m1+m2)/2；
（5）重复步骤（2）～（4），直到阈值变化小于设定值。
示例：迭代算法分割带噪声的指纹图像
void threshold_basic(Mat& src, Mat& dst) {
int totalGray = sum(src)[0];
int totalPixels = src.rows * src.cols;
int T_pre = totalGray / totalPixels; double T = 0; double dT = T_pre - T;
while (abs(dT) > 1) {
int G1 = 0; int G2 = 0;
int count1 = 0; int count2 = 0;
for (int i = 0; i < src.rows; i++) {
for (int j = 0; j < src.cols; j++) {
int gray = src.at<uchar>(i, j);
if (gray > T_pre) {
count1++;
G1 += gray;
}
else {
count2++;
G2 += gray;
}
}
}
T = (G1 * 1.0 / count1 + G2 * 1.0 / count2) / 2;
dT = T_pre - T;
T_pre = T;
}
threshold(src, dst, T, 255, THRESH_BINARY);
}
3 使用Otsu方法的最优全局阈值处理
阈值处理本质上是对像素进行分类的统计决策问题，其目的是在把像素分配给两组或者多组（也称分类）的过程中，使引入的平均误差最小。OTSU 方法又称大津算法，使用最大化类间方差（intra-class variance）作为评价准则，基于对图像直方图的计算，可以给出类间最优分离的最优阈值。任取一个灰度值 T，可以将图像分割为两类C1和C2，C1和C2的像素数的占比分别为 P1、P2，C1和C2的灰度值均值分别为 m1、m2，整个图像的平均灰度为 m，定义类间方差为：
为了评估阈值的有效性，我们使用归一化的、无量纲的可分离性测度：
目标是求最大化类间方差的阈值k。
示例：使用Otsu方法的最优全局阈值处理
下面是聚合细胞的光学显微图像。目的是从背景中分割出分子。由于直方图没有明显的波谷，并且背景和目标之间的灰度差很小，上面的基于全局的迭代算法没有实现期望的分割。使用基本算法算出的阈值是169，Otsu方法计算得到的阈值是182，可分离性测度为0.467。
// opencv函数调用
threshold(src, dst, 0, 255, THRESH_OTSU);
// Otsu 算法实现
double threshold_otsu(Mat& src, Mat& dst, int& k) {
// 生成直方图 int histsize = 256; float range[] = { 0, 256 }; const float* ranges = { range }; Mat hist; calcHist(&src, 1, 0, Mat(), hist, 1, &histsize, &ranges); drawHist(hist, dst); int totalGray = sum(src)[0]; int totalPixels = src.rows * src.cols;
// 计算全局灰度均值 double m = totalGray * 1.0 / totalPixels; double C_std = 0; double G_std; for (int T = 100; T < 200; T++) { int numC1 = 0; int sumC1 = 0; int numC2 = 0; int sumC2 = 0; G_std = 0; double k_sum = 0; double k_count=0; for (int i = 1; i < 256; i++) { if (i < T) { numC1 += hist.at<float>(i, 0); sumC1 += hist.at<float>(i, 0) * i; } else { numC2 += hist.at<float>(i, 0); sumC2 += hist.at<float>(i, 0) * i; } G_std += pow(i - m, 2) * hist.at<float>(i, 0); } if (numC1 == 0) continue; double P1 = numC1 * 1.0 / totalPixels; double m1 = sumC1 * 1.0 / numC1; double m2 = sumC2 * 1.0 / numC2; double P2 = 1 - P1;
// 计算类间方差 double tmp = P1 * pow(m1 - m, 2) + P2 * pow(m2 - m, 2); if (tmp >= C_std) { if (tmp > C_std) { k_sum = T; k_count = 1; C_std = tmp; k = T; } else { //如果极大值不唯一，取极大值的平均 k_sum += T; k_count++; k = k_sum / k_count * 1.0; } } } threshold(src, dst, k, 255, THRESH_BINARY); return C_std / G_std * totalPixels; // 返回可分离性测度 }
4 使用图像平滑改进全局阈值处理
噪声会使得简单的阈值处理问题变得不可求解。无法从源头降低噪声，并且阈值处理是首选分割方法时，增强性能的一种技术通常是在阈值处理之前平滑图像。
下图(a)是带噪声的图像，图(b)显示了其直方图，图(c)时Otsu方法对图像进行阈值处理后的结果。白色区域中的每个黑点和黑色区域中的每个白点都是阈值处理的误差，所以分割的很不成功。图(d)是用一个5×5的平均核平滑带噪声图像后的结果，图(e)是其直方图。由于平滑处理，直方图形状的改进很明显，因此阈值处理后的结果也很完美，如图(f)。但是，分割后的目标和背景之间的边界有轻微的失真，这是由于平滑图像造成的边界模糊。对图像平滑的次数越多，分割后的结果的边界误差越大。
5 使用边缘改进全局阈值处理
上个例子体现了图像平滑对阈值处理的好处，但是当前景区域相对于背景的尺寸大幅缩小时，利用5×5的平均核平滑后，直方图没有清晰的波谷，分割还是会失败。也就是说，前景区域小到无法与噪声引起的灰度扩展相比，对直方图的贡献很小。
如果只利用位于或接近目标和背景间的边缘的像素，那么得到的直方图将有几个高度近似相同的波峰。此外任何位于目标上的像素的概率将近似等于其位于背景上的概率，从而改进了直方图模式的对称性。如何选取这些像素呢，我们可以通过计算其梯度或拉普拉斯来确定。例如，在边缘的过渡处，拉普拉斯的平均值为零，因此根据拉普拉斯准则选取的像素所形成的直方图的波谷应该很少。这一性质往往会使得直方图产生较深的波谷。在实践中，通常使用梯度图像或拉普拉斯图像来得到类似的结果，后者更受青睐，因为它在计算上更有吸引力，而且也是使用各向同性边缘检测子创建的。前述讨论可小结为如下算法，其中f(x,y)是输入图像
1. 使用任意梯度算子，计算f(x,y)的梯度幅度或拉普拉斯的绝对值
2. 规定一个阈值，对上一步的图像进行阈值处理，得到强边缘像素作为模板
3. 只用f(x,y)中对应于g(x,y)中1值像素位置的那些像素计算直方图
4. 用上一步中的直方图全局地分割f(x,y)，例如使用Otsu
示例：
Mat src = imread("./12.tif", 0);
// 计算梯度
Mat gradX, gradY;
Sobel(src, gradX, CV_32F, 1, 0, 3);
Sobel(src, gradY, CV_32F, 0, 1, 3);
convertScaleAbs(gradX, gradX);
convertScaleAbs(gradY, gradY);
Mat grad = gradX + gradY;
double minValue, maxValue;
Point minIdx, maxIdx;
minMaxLoc(grad, &minValue, &maxValue, &minIdx, &maxIdx);
// 阈值处理获得边缘
double T = maxValue * 0.75;
Mat mask;
threshold(grad, mask, T, 255, THRESH_BINARY);
Mat edgeImg = Mat::zeros(src.size(), CV_8UC1);
src.copyTo(edgeImg, mask);
// 对于强边缘计算直方图
int histsize = 256;
float range[] = { 0, 256 };
const float* ranges = { range };
Mat hist;
calcHist(&edgeImg, 1, 0, Mat(), hist, 1, &histsize, &ranges);
Mat dst = Mat::zeros(256, 256, CV_8UC3);
for (int i = 0; i < 255; i++) {
int h = (int)hist.at<float>(i, 0);
if (h > 0)
rectangle(dst, Rect(i, 256 - h, 1, h), Scalar(255,255,255), -1);
}
// 观察直方图得阈值为135左右 （用otsu找不到那个阈值==！）
rectangle(dst, Rect(135, 0, 1, 256), Scalar(0,255,0), -1);
Mat binImg;
threshold(src, binImg, 135, 255, THRESH_BINARY);
示例：下图是酵母细胞的一幅8比特图像，我们希望利用全局阈值处理来得到图像中与亮点对应的区域。
Otsu方法算出的阈值是42，可分离性测度是0.636。对强边缘计算阈值为144，可分离性测度是0.98193。
void threshold_grad(Mat& src, Mat& dst, double percent) {
Mat grad;
Laplacian(src, grad, CV_64F, 3);
convertScaleAbs(grad, grad);
double minValue, maxValue;
Point minIdx, maxIdx;
minMaxLoc(grad, &minValue, &maxValue, &minIdx, &maxIdx);
double T = maxValue * percent;
Mat mask;
threshold(grad, mask, T, 255, THRESH_BINARY);
Mat edgeImg = Mat::zeros(src.size(), CV_8UC1);
src.copyTo(edgeImg, mask);
Mat tmp; int k;
double n = threshold_otsu(edgeImg, tmp, k);
threshold(src, dst, 130, 255, THRESH_BINARY);
}
int main() {
Mat src = imread("./13.tif", 0);
Mat binImg; int k;
// 1. Otsu
double n = threshold_otsu(src, binImg, k); // 得到k=42， n=0.635996
// 2. 梯度+Otsu
threshold_grad(src, binImg, 0.7);
}
6 多阈值处理
迄今为止，我们关注的是用单个阈值对图像进行分割。Otsu方法可以扩展到任意数量的阈值，因为基于这种方法的可分离性测度也可拓展到任意数量的类。但是随着类数的增加，它开始失去意义，因为我们只处理一个变量（灰度）。要求两个以上阈值的应用，通常使用更多的灰度值来求解。因此，下面讨论双阈值处理。对于由三个灰度区间组成的三个类（三个类由两个阈值分隔），类间方差由下式给出：
其中，m1，m2，m3分别是三个类别的平均灰度值，P1，P2，P3分别是三类像素占总像素的比，mG为整张图的平均灰度。使类间方差 最大化的灰度值 k1和k2就是最优阈值。阈值处理后的图像由下式给出：
可分离性测度可直接拓展到多个阈值：
示例：一幅冰山的图像，分割为暗背景、冰山的明亮区域和阴影区域。
从直方图可以看出，需要两个阈值。计算得到，两个阈值分别为T1=80，T2=177，可分离性测度是0.954。分割得这样好的主要原因是，能够找到直方图中由合适宽度和深度的波谷分隔的3个不同模式。
double threshold_otsu1(Mat& src, Mat& dst, int& T1_final, int& T2_final) {
int histsize = 256;
float range[] = { 0, 256 };
const float* ranges = { range };
Mat hist;
calcHist(&src, 1, 0, Mat(), hist, 1, &histsize, &ranges);
int totalGray = sum(src)[0];
int totalPixels = src.rows * src.cols;
double m = totalGray * 1.0 / totalPixels;
double C_std = 0; double G_std;
for (int T1 = 1; T1 < 255; T1++) {
for (int T2 = 2; T2 < 255; T2++) {
int numC1 = 0; int sumC1 = 0;
int numC2 = 0; int sumC2 = 0;
int numC3 = 0; int sumC3 = 0;
G_std = 0; double k1_sum = 0; double k2_sum = 0; double k_count = 0;
for (int i = 0; i < 256; i++) {
if (i <= T1) {
numC1 += hist.at<float>(i, 0);
sumC1 += hist.at<float>(i, 0) * i;
}
else if (i <= T2) {
numC2 += hist.at<float>(i, 0);
sumC2 += hist.at<float>(i, 0) * i;
}
else {
numC3 += hist.at<float>(i, 0);
sumC3 += hist.at<float>(i, 0) * i;
}
G_std += pow(i - m, 2) * hist.at<float>(i, 0);
}
if (numC1 == 0 || numC2 == 0 || numC3 == 0)
continue;
double P1 = numC1 * 1.0 / totalPixels;
double P2 = numC2 * 1.0 / totalPixels;
double m1 = sumC1 * 1.0 / numC1;
double m2 = sumC2 * 1.0 / numC2;
double m3 = sumC3 * 1.0 / numC3;
double P3 = 1 - P1 - P2;
double tmp = P1 * pow(m1 - m, 2) + P2 * pow(m2 - m, 2) + P3 * pow(m3 - m, 2);
if (tmp >= C_std) {
if (tmp > C_std) {
k1_sum = T1;
k2_sum = T2;
k_count = 1;
C_std = tmp;
T1_final = T1;
T2_final = T2;
}
else {
k1_sum += T1;
k2_sum += T2;
k_count++;
T1_final = (int)(k1_sum / k_count * 1.0);
T2_final = (int)(k2_sum / k_count * 1.0);
}
}
}
}
dst = Mat(src.size(), CV_8UC3);
vector<Mat> bgr;
split(dst, bgr);
for (int i = 0; i < src.rows; i++) {
for (int j = 0; j < src.cols; j++) {
int gray = src.at<uchar>(i, j);
if (gray <= T1_final)
bgr[0].at<uchar>(i, j) = 255;
else if (gray <= T2_final)
bgr[1].at<uchar>(i, j) = 255;
else
bgr[2].at<uchar>(i, j) = 255;
}
}
merge(bgr, dst);
return C_std / G_std * totalPixels;
}
7 可变阈值处理
前面的全局阈值算法对于某些光照不均的图像会束手无策，接下来将介绍局部阈值算法。其基本方法是在图像中的每个点，根据其邻域的一条或多条规定性质计算阈值。标准差和均值对于求局部阈值是对比度和平均灰度的描述子，在局部阈值处理中非常有用。令mxy和σxy是图像中坐标(x,y)为中心的邻域Sxy所包含的像素集的标准差和均值。下面是基于局部性质的可变阈值的通用公式：
分割后的图像计算为
式中，Q是一个谓词逻辑，它是用邻域Sxy中的像素算出的谓词逻辑，可给可变阈值处理分配更大的权重：
示例：使用局部阈值处理图像
下面是一幅被斑点灰度模式遮蔽的手写文本图像。这种灰度遮蔽在用点照明得到的图像中很常见。Otsu全局处理方法的效果较差，而利用n=20，c=0.5的局部阈值效果很好。
void threshold_adaptive(Mat& src, Mat& dst, int ksize, double c) {
Mat mean;
blur(src, mean, Size(ksize, ksize));
dst = Mat::zeros(src.size(), CV_8UC1);
for (int i = 0; i < src.rows; i++) {
for (int j = 0; j < src.cols; j++) {
if (src.at<uchar>(i, j) > c * mean.at<uchar>(i, j))
dst.at<uchar>(i, j) = 255;
}
}
}
...
threshold_adaptive(src, binImg, 20, 0.5);//调用
adaptiveThreshold(src, dst, 255, ADAPTIVE_THRESH_GAUSSIAN_C, THRESH_BINARY, 5, 3);//使用opencv函数
OpenCV提供了函数cv::adaptiveThreshold来实现自适应阈值处理。
void adaptiveThreshold( InputArray src, OutputArray dst, double maxValue, int adaptiveMethod, int thresholdType, int blockSize, double C );
scr：输入图像，nparray 二维数组，必须是 8-bit 单通道灰度图像！
dst：输出图像，大小和格式与 scr 相同
maxValue：为满足条件的像素指定的非零值，具体用法见 thresholdType 说明
adaptiveMethod：自适应方法选择，其中ADAPTIVE_THRESH_MEAN_C阈值是邻域的均值；ADAPTIVE_THRESH_GAUSSIAN_C阈值是邻域的高斯核加权均值。
thresholdType：阈值处理方法，THRESH_BINARY：大于阈值时置 maxValue，否则置 0；THRESH_BINARY_INV：大于阈值时置 0，否则置 maxValue
blockSize：像素邻域的尺寸，用于计算邻域的阈值，通常取 3，5，7
C：偏移量，从邻域均值中减去该常数
参考
1. 冈萨雷斯《数字图像处理（第四版）》Chapter 10（所有图片可在链接中下载）
