二维Otsu算法的原理与实现

1.简介：

       一维Otsu算法有计算简洁、稳定、自适应强等优点，被广泛用于图像分割中。但一维Otsu算法没有考虑图像像素点之间的关系，当图像中有噪声时，会导致分割的效果不理想。因此，刘健庄等人在1993年提出了二维的Otsu算法，提升了算法的抗噪声能力。

2.算法思想：

       同时考虑像素的灰度值分布和它们邻域像素的平均灰度值分布，因此形成的阈值是一个二维矢量，最佳的阈值在一个二维的测度准则下确定最大值时得到。

3.算法过程：

（1）设图像I(x,y),的灰度级为L级，那么图像的邻域平均灰度也分为L级。
（2）设f(x,y)为像素点(x,y)的灰度值，g(x,y)为像素点(x,y)为中心的K*K的像
        素点集合的灰度平均值。令f(x,y)=i,g(x,y)=j，然后就形成了一个二元组(i,j)。
（3）设二元组（i,j）出现的次数为fij，然后求出二元组对应的概率密度Pij，
        Pij=fij/N, i,j=1,2，…，L，其中N为图像像素点总数。
（4）任意选取一个阈值向量(s,t)选取的阈值向量将图像的二维直方图划分成4个
        区域，B、C区域代表图像的前景和背景，A、D区域代表噪声点。

（5）设C、B两个区域对应的概率分别为w1，w2，对应的均值矢量为u1，u2。整个图
        片所对应的均值矢量为uT。

4.代码实现（opencv3）：

    

#include "opencv.hpp"#include "imgproc.hpp"#include "highgui.hpp"#include "iostream"#include "core.hpp"using namespace cv;using namespace std;int Otsu2D(Mat srcimage);         //二维Otsu算法int main(){Mat srcimage, grayimage, dstimage;srcimage = imread("lena.jpg");namedWindow("原图", 0);imshow("原图", srcimage);                                            //显示原图cvtColor(srcimage, grayimage, COLOR_RGB2GRAY);                         //得到灰度图double time0 = static_cast<double>(getTickCount());                    //记录程序开始时间int thresholdValue = Otsu2D(grayimage);                              //调用二维Otsu函数time0 = ((double)getTickCount() - time0) / getTickFrequency();cout << "算法运行时间为：" << time0 << endl;cout << "Otsu阈值为：" << thresholdValue << endl;threshold(grayimage, dstimage, thresholdValue, 255, THRESH_BINARY);    //将得到的阈值传入函数，得到分割效果图namedWindow("Otsu算法结果", 0);imshow("Otsu算法结果", dstimage);waitKey();    return 0;}

int Otsu2D(Mat srcimage){double Histogram[256][256];        //建立二维灰度直方图double TrMax = 0.0;                //用于存储矩阵的迹（矩阵对角线之和）int height = srcimage.rows;        //矩阵的行数int width = srcimage.cols;         //矩阵的列数int N = height*width;              //像素的总数int T;                             //最终阈值uchar *data = srcimage.data;for (int i = 0; i < 256; i++){for (int j = 0; j < 256; j++){Histogram[i][j] = 0;      //初始化变量}}for (int i = 0; i < height; i++){for (int j = 0; j < width; j++){int Data1 = data[i*srcimage.step + j];         //获取当前灰度值int Data2 = 0;                           //用于存放灰度的平均值for (int m = i - 1; m <= i + 1; m++){for (int n = j - 1; n <= j + 1; n++){if ((m >= 0) && (m < height) && (n >= 0) && (n < width))Data2 += data[m*srcimage.step + n];//邻域灰度值总和}}Data2 = Data2 / 9;Histogram[Data1][Data2]++;                  //记录（i,j）的数量}}for (int i = 0; i < 256; i++)for (int j = 0; j < 256; j++)Histogram[i][j] /= N;     //归一化的每一个二元组的概率分布double Fgi = 0.0;    //前景区域均值向量i分量double Fgj = 0.0;    //前景区域均值向量j分量double Bgi = 0.0;    //背景区域均值向量i分量double Bgj = 0.0;    //背景区域均值向量j分量double Pai = 0.0;    //全局均值向量i分量 panorama(全景)double Paj = 0.0;    //全局均值向量j分量double w0 = 0.0;     //前景区域联合概率密度double w1 = 0.0;     //背景区域联合概率密度double num1 = 0.0;   //遍历过程中前景区i分量的值double num2 = 0.0;   //遍历过程中前景区j分量的值double num3 = 0.0;   //遍历过程中背景区i分量的值double num4 = 0.0;   //遍历过程中背景区j分量的值int Threshold_s = 0; //阈值sint Threshold_t = 0; //阈值tdouble temp = 0.0;   //存储矩阵迹的最大值for(int i=0;i<256;i++){for (int j = 0; j < 256; j++){Pai += i*Histogram[i][j];   //全局均值向量i分量计算Paj += j*Histogram[i][j];   //全局均值向量j分量计算}}for (int i = 0; i < 256; i++){for (int j = 0; j < 256; j++){w0 += Histogram[i][j];        //前景的概率num1 += i*Histogram[i][j];    //遍历过程中前景区i分量的值num2 += j*Histogram[i][j];    //遍历过程中前景区j分量的值w1 = 1 - w0;                  //背景的概率num3 = Pai - num1;            //遍历过程中背景区i分量的值num4 = Paj - num2;            //遍历过程中背景区j分量的值Fgi = num1 / w0;              Fgj = num2 / w1;Bgi = num3 / w0;Bgj = num4 / w1;TrMax = ((w0*Pai - num1)*(w0*Pai - num1) + (w0*Paj - num2)*(w0*Paj - num2)) / (w0*w1);if (TrMax > temp){temp = TrMax;Threshold_s = i;Threshold_t = j;}}}cout << Threshold_s << " " << Threshold_t << endl;T = Threshold_s;return T;}

5.运行结果（得到2个阈值，取其中一个即可）：