6.1、方法概述

阈值分割的核心就是如何选取阈值， 选取正确的阈值是分割成功的关键。

1、全局阈值分割

全局阈值分割指的是将灰度值大于thresh（阈值）的像素设为白色，小于或者等于 thresh的像素设为黑色； 或者反过来， 将大于thresh的像素设为黑色， 小于或者等于thresh 的像素设为白色， 两者的区别只是呈现形式不同。

image.png

double threshold( InputArray src, OutputArray dst,double thresh, double maxval, int type );

image.png

需要注意的是，当类型为THRESH_OTSU或THRESH_TRIANGLE时，输入参数src只支持uchar类型， 这时thresh也是作为输出参数的， 即通过Otsu和TRIANGLE算法自动计算出来。

2、局部阈值分割

局部阈值分割的核心也是计算阈值矩阵，比较常用的是后面提到的自适应阈值算法（又称移动平均值算法） ， 是一种简单但是高效的局部阈值算法，其核心思想就是把每一个像素的邻域的“平均值”作为该位置的阈值。

6.2、直方图技术法

一幅含有一个与背景呈现明显对比的物体的图像具有包含双峰的直方图，两个峰值对应于物体内部和外部较多数目的点，两个峰值之间的波谷对应于物体边缘附近相对较少数目的点。

直方图技术法就是首先找到这两个峰值，然后取两个峰值之 间的波谷位置对应的灰度值，就是所要的阈值。

一 种常用的方式是先对直方图进行高斯平滑处理，逐渐增大高斯滤波器的标准差，直到能从平滑后的直方图中得到两个唯一的波峰和它们之间唯一的最小值。但这种方式需要手动调节，下面介绍一种规则自动选取波峰和波谷的方式。

假设输入图像为I， 高为H、 宽为W，histogramI代表其对应的灰度直方图， histogramI （k） 代表灰度值等于k的像素点个数， 其中0≤k≤255。

• 第一步： 找到灰度直方图的第一个峰值，并找到其对应的灰度值。显然，灰度直方图的最大值就是第一个峰值且对应的灰度值用firstPeak表示。

• 第二步： 找到直方图的第二个峰值，并找到其对应的灰度值。第二个峰值不一定是直方图的第二大值，因为它很有可能出现在第一个峰值的附近。

image.png

• 第三步： 找到这两个峰值之间的波谷，如果出现两个或者多个波谷，则取左侧的波 谷即可，其对应的灰度值即为阈值。

int threshTwoPeaks(const Mat &image, Mat &thresh_out)
{
    //计算灰度直方图
    Mat histogram = calGrayHist(image);
    //找到灰度直方图最大峰值对应的灰度值
    Point firstPeakLoc;
    minMaxLoc(histogram, NULL, NULL, NULL, &firstPeakLoc);
    int firstPeak = firstPeakLoc.x;
    //寻找灰度直方图第二个峰值对应的灰度值
    Mat measureDists = Mat::zeros(Size(256, 1), CV_32FC1);
    for (int k = 0; k < 256; k++)
    {
        int hist_k = histogram.at<int>(0, k);
        measureDists.at<float>(0, k) = pow(float(k - firstPeak), 2)*hist_k;
    }
    Point secondPeakLoc;
    minMaxLoc(measureDists, NULL, NULL, NULL, &secondPeakLoc);
    int secondPeak = secondPeakLoc.x;

    //找到两个之间的最小值对应的灰度值，作为阈值
    Point threshLoc;
    int thresh = 0;
    if (firstPeak < secondPeak)//第一个峰值在第二个峰值的左侧
    {
        minMaxLoc(histogram.colRange(firstPeak, secondPeak), NULL, NULL, NULL, &threshLoc);
        thresh = firstPeak + threshLoc.x + 1;
    }
    else//第一个峰值在第二个峰值的右侧
    {
        minMaxLoc(histogram.colRange(secondPeak, firstPeak), NULL, NULL, NULL, &threshLoc);
        thresh = secondPeak + threshLoc.x + 1;
    }
    //阈值分割
    threshold(image, thresh_out, thresh, 255, THRESH_BINARY);
    return thresh;
}

6.3、熵方法

image.png

利用熵计算阈值的步骤如下:

• 第一步： 计算I 的累加概率直方图，又称零阶累积矩，记为

image.png

• 第二步： 计算各个灰度级的熵，记为

image.png

• 第三步： 计算使f （t） =f1（t） +f2（t） 最大化的t值， 该值即为得到的阈值， 即thresh=argtmax（f （t） ）， 其中

image.png

6.4、Otsu阈值处理

在对图像进行阈值分割时，所选取的分割阈值应使前景区域的平均灰度、背景区域 的平均灰度与整幅图像的平均灰度之间的差异最大， 这种差异用区域的方差来表示。 Otsu[2]提出了最大方差法， 该算法是在判别分析最小二乘法原理的基础上推导得出的， 计算过程简单， 是一种常用的阈值分割的稳定算法。

• 第一步： 计算灰度直方图的零阶累积矩（或称累加直方图） 。

image.png

• 第二步： 计算灰度直方图的一阶累积矩。

image.png

• 第三步： 计算图像I 总体的灰度平均值mean， 其实就是k=255时的一阶累积距， 即

image.png

• 第四步： 计算每一个灰度级作为阈值时， 前景区域的平均灰度、 背景区域的平均灰 度与整幅图像的平均灰度的方差。 对方差的衡量采用以下度量：

image.png

• 第五步： 找到上述最大的σ2（k）， 然后对应的k即为Otsu自动选取的阈值， 即

image.png

int otsu(const Mat &image, Mat &OtsuThreshImage)
{
    //计算灰度直方图
    Mat histogram = calGrayHist(image);
    //归一化灰度直方图
    Mat normHist;
    histogram.convertTo(normHist, CV_32FC1, 1.0 / (image.rows*image.cols), 0.0);
    //计算累加直方图（零阶累积距）和一阶累积距
    Mat zeroCumuMoment = Mat::zeros(Size(256, 1), CV_32FC1);
    Mat oneCumuMoment = Mat::zeros(Size(256, 1), CV_32FC1);
    for (int i = 0; i < 256; i++)
    {
        if (i == 0)
        {
            zeroCumuMoment.at<float>(0, i) = normHist.at<float>(0, i);
            oneCumuMoment.at<float>(0, i) = normHist.at<float>(0, i);
        }
        else
        {
            zeroCumuMoment.at<float>(0, i) = normHist.at<float>(0, i) + zeroCumuMoment.at<float>(0, i - 1);
            oneCumuMoment.at<float>(0,i)= normHist.at<float>(0, i) + oneCumuMoment.at<float>(0, i - 1);
        }
    }
    //计算类间方差
    Mat variance = Mat::zeros(Size(256, 1), CV_32FC1);
    //总平均值
    float mean = oneCumuMoment.at<float>(0, 255);
    for (int j = 0; j < 256; j++)
    {
        if (zeroCumuMoment.at<float>(0, j) == 0 || zeroCumuMoment.at<float>(0, j) == 1)
        {//分母不能为零
            variance.at<float>(0, j) = 0;
        }
        else
        {
            variance.at < float>(0, j) = pow(mean*zeroCumuMoment.at<float>(0, j) - oneCumuMoment.at<float>(0, j), 2) / zeroCumuMoment.at<float>(0, j)*(1.0 - zeroCumuMoment.at<float>(0, j));
        }
    }
    //找到最值作为阈值
    Point maxLoc;
    minMaxLoc(variance, NULL, NULL, NULL, &maxLoc);
    int thresh = maxLoc.x;
    //阈值处理
    threshold(image, OtsuThreshImage, thresh, 255, THRESH_BINARY);
    return thresh;
}

6.5、自适应阈值

在不均匀照明或者灰度值分布不均的情况下，如果使用全局阈值分割， 那么得到的分割效果往往会很不理想。那么想到的策略是针对每一个位置的灰度值 设置一个对应的阈值， 而该位置阈值的设置也和其邻域有必然的关系。

在对图像进行平滑处理时，均值平滑、高斯平滑、中值平滑用不同规则计算出以当前像素为中心的邻域内的灰度“平均值”， 所以可以使用平滑处理后的输出结果作为每个 像素设置阈值的参考值，如用均值滤波后的结果乘以某个比例系数作为最后的阈值矩阵。

平滑算子的宽度必须大于被识别物体的宽度，平滑算子的尺寸越大，平滑后的结果越能更好地作为每个像素的阈值的参考，当然也不能无限大。

• 第一步： 对图像进行平滑处理， 平滑结果记为fsmooth（I） ，其中fsmooth可以代表 均值平滑、高斯平滑、中值平滑。
• 第二步： 自适应阈值矩阵Thresh=（1-ratio） *fsmooth（I） ，一般令ratio=0.15。
• 第三步： 利用局部阈值分割的规则进行阈值分割。

Mat adaptiveThresh(Mat I, int radius, float radio, METHOD method)
{
    //step1:对图像矩阵进行平滑处理
    Mat smooth;
    switch (method)
    {
    case MEAN:
        boxFilter(I, smooth, CV_32FC1, Size(2 * radius + 1, 2 * radius + 1));
        break;
    case GAUSS:
        GaussianBlur(I, smooth, Size(2 * radius + 1, 2 * radius + 1), 0, 0);
        break;
    case MEDIAN:
        medianBlur(I, smooth, 2 * radius + 1);
        break;
    default:
        break;
    }
    //step2:平滑结果乘于比例系数，然后图像矩阵与其做差
    I.convertTo(I, CV_32FC1);
    smooth.convertTo(smooth, CV_32FC1);
    Mat diff = I - (1.0 - radio)*smooth;
    //step3:阈值处理，当≥0，输出值为255，反之，输出值为0
    Mat out = Mat::zeros(diff.size(), CV_8UC1);
    for (int r = 0; r < I.rows; r++)
    {
        for (int c = 0; c < I.cols; c++)
        {
            if (diff.at<float>(r, c) > 0)
            {
                out.at<uchar>(r, c) = 255;
            }
        }
    }
    return out;
}

就可以理解OpenCV提供的自适应阈值函数：

void adaptiveThreshold( InputArray src, OutputArray dst,
                                     double maxValue, int adaptiveMethod,
                                     int thresholdType, int blockSize, double C );

6.6、二值图的逻辑运算

OpenCV提供的两个函数bitwise_and和bitwise_or分别实现了两 个矩阵之间的与运算和或运算，它们本质上完成的是两个矩阵对应位置数值的逻辑运算。