版权声明：本文为博主原创文章，未经博主允许不得转载。 https://blog.csdn.net/kksc1099054857/article/details/78321776
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;
-
}
- <code class="language-cpp">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; //阈值s
- int Threshold_t = 0; //阈值t
- double 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;
- }
- </code>
-
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;
//阈值s
-
int Threshold_t =
0;
//阈值t
-
double 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个阈值，取其中一个即可）：
6.修改后的代码
-
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 S[
256];
//统计前i行概率的数组
-
double N1[
256];
//统计遍历过程中前景区i分量的值
-
double N2[
256];
//统计遍历过程中前景区j分量的值
-
S[
0] =
0;
-
N1[
0] =
0;
-
N2[
0] =
0;
-
for (
int i =
1; i <
256; i++)
-
{
-
double x =
0,n1=
0,n2=
0;
-
for (
int j =
0; j <
256; j++)
-
{
-
x += Histogram[i -
1][j];
-
n1 += ((i
-1)*Histogram[i -
1][j]);
//遍历过程中前景区i分量的值
-
n2 += (j*Histogram[i -
1][j]);
//遍历过程中前景区j分量的值
-
}
-
S[i] = x + S[i -
1];
-
N1[i] = n1 + N1[i -
1];
-
N2[i] = n2 + N2[i -
1];
-
}
-
-
double Pai =
0.0;
//全局均值向量i分量 panorama(全景)
-
double Paj =
0.0;
//全局均值向量j分量
-
int Threshold_s =
0;
//阈值s
-
int Threshold_t =
0;
//阈值t
-
int M =
0;
//中间变量
-
double temp =
0.0;
//存储矩阵迹的最大值
-
double num3 =
0.0;
//遍历过程中背景区i分量的值
-
double num4 =
0.0;
//遍历过程中背景区j分量的值
-
double Fgi =
0.0;
//前景区域均值向量i分量
-
double Fgj =
0.0;
//前景区域均值向量j分量
-
double Bgi =
0.0;
//背景区域均值向量i分量
-
double Bgj =
0.0;
//背景区域均值向量j分量
-
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++)
-
{
-
double w0 =
0.0;
//前景区域联合概率密度
-
double w1 =
0.0;
//背景区域联合概率密度
-
double num1 =
0.0;
//遍历过程中前景区i分量的值 与w0一样做相关处理
-
double num2 =
0.0;
//遍历过程中前景区j分量的值
-
-
if (i >=
1)
-
{
-
w0 += S[i -
1];
-
num1 += N1[i -
1];
-
num2 += N2[i -
1];
-
}
-
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;
-
}
-
}
-
cout << Threshold_s <<
endl;
-
T = Threshold_s;
-
return T;
-
}
扫一扫
6019
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
