重庆分公司
重庆分公司，新征程启航
为企业提供网站建设、域名注册、服务器等服务
C++手敲基于梯度图和像素数量数组的OTSU阈值分割-创新互联
一、OTSU算法原理
文章标题：C++手敲基于梯度图和像素数量数组的OTSU阈值分割-创新互联
本文网址：http://cqcxhl.cn/article/dioeed.html
➢OTSU法（大类间方差法，有时也称之为大津算法）
➢ 使用聚类的思想，把图像的灰度数按灰度级分成2个部分， 使得两个部分之间的灰度值差异大，每个部分之间的灰 度差异最小
➢ 通过方差的计算来寻找一个合适的灰度级别来划分。
➢ 可以在二值化的时候采用OTSU算法来自动选取阈值进行 二值化。
➢ OTSU算法被认为是图像分割中阈值选取的最佳算法，计 算简单，不受图像亮度和对比度的影响。
➢ 因此,使类间方差大的分割意味着错分概率最小。
二、OTSU算法步骤：全局阈值T可以按如下计算：
➢选择一个初时估计值T（一般为图像的平均灰度值）
➢使用T分割图像，产生两组像素：G1包括灰度级大于T的像素，G2包括灰度级小于等于T的像素
➢计算G1 中像素的平均值并赋值给μ1，计算G2 中像素 的平均值并赋值给μ2
➢计算一个新的阈值：
➢重复步骤 2 ～ 4，一直到两次连续的T之间的差小于 预先给定的上界T∞
三、代码实现当我们用边缘梯度算子（如Roberts、Soberl、Prewitt等）处理原始图像拿到梯度图之后，可以创建一个1*256维的数组来存放梯度图中每个像素的个数，其具体的实现方式如下：
int pixel_num[256] = { 0 }; //数组记得初始化为0，即未统计数量之前各个像素的个数都是0
for (int row = 0; row< img.rows; row++) {
for (int col = 0; col< img.cols; col++) {
pixel_num01[Sobel_City.at(row, col)] += 1; //遍历到的像素值作为索引，次数+1
}
} 以下是OTSU的实现方法：
①由传入的梯度图拿到图像的维度，创建一个空白图像用于保留运算得到的结果，并且拿到梯度图像的总像素个数。
Mat OTSU = Mat::zeros(Size(img.cols, img.rows), img.type());
int N = img.rows * img.cols;②由像素数量数组，得到每一个像素对应的概率：
double pro[256] = { 0 }; // 定义概率数组
for (int i = 0; i< 256; i++) {
pro[i] = arr[i] / N; // 得到每个像素值的概率
}注意创建的概率数组是浮点型，定义为整形的话会出现全部像素的概率都为零的情况。
③依次创建OTSU计算的过程量，包括两类的概率、均值，还有阈值T以及方差v等等；先计算得到两类的均值以及概率，根据公式v = w0 * w1 * pow(u1 - u0, 2);就可以得到此次遍历的像素作为阈值得到的方差。
double delta = 0, w0 = 0, w1 = 0, u0 = 0, u1 = 0, v = 0; // 变量初始化为0
int T = 0, thresh = 0;
while (T<= 255) {
for (int i = 0; i<= T; i++) {
w0 += pro[i]; // C0的概率
}
w1 = 1 - w0; // C1的概率
for (int i = 0; i<= 255; i++) {
if (i<= T) {
u0 += pro[i] * i; // C0的平均值
}
else {
u1 += pro[i] * i; // C1的平均值
}
v = w0 * w1 * pow(u1 - u0, 2);
if (v >delta) {
delta = v;
thresh = T;
}
T += 1;
}
}遍历循环的过程中要对方差进行判断，如果此次计算的方差要比之前的大方差delta要大，那么对delta进行更新，同时也用thresh记录下历史大方差对应的像素值，作为分割的阈值T。这样一来，遍历的结果就可以得到我们想要的，能使两类差异大的阈值T。
④通过对与梯度图同维大小的空白图片遍历，给每一个像素进行二值化。若遍历到的像素值大于等于T，则赋值为255，否则为0。
for (int row = 0; row< img.rows; row++) {
for (int col = 0; col< img.cols; col++) {
if (img.at(row, col) >thresh) { //遍历判断
OTSU.at(row, col) = 255;
}
}
} ⑤显示OTSU图像，得到阈值分割结果：
imshow("OTSU阈值分割", OTSU);
waitKey(0);
destroyAllWindows();
四、代码汇总
#include#include#includeusing namespace cv;
using namespace std;
void Sobel(Mat img);
void OTSU(Mat& img, int arr[]);
int main() {
Mat Gray = imread("C:\\Users\\Czhannb\\Desktop\\gray.png", IMREAD_GRAYSCALE);
if (Gray.empty()) {
cout<< "读取图片错误！"<< endl;
}
else {
imshow("未动工之前：", Gray);
}
Sobel(Gray);
return 0;
}
void OTSU(Mat& img, int arr[]) { // 需要输入待处理的图像 和 直方图像素个数数组
int N = img.rows * img.cols; // 统计输入图像的像素个数N
double pro[256] = { 0 }; // 定义概率数组
for (int i = 0; i< 256; i++) {
pro[i] = arr[i] / N; // 得到每个像素值的概率
}
double delta = 0, w0 = 0, w1 = 0, u0 = 0, u1 = 0, v = 0; // 变量初始化为0
int T = 0, thresh = 0;
while (T<= 255) {
for (int i = 0; i<= T; i++) {
w0 += pro[i]; // C0的概率
}
w1 = 1 - w0; // C1的概率
for (int i = 0; i<= 255; i++) {
if (i<= T) {
u0 += pro[i] * i; // C0的平均值
}
else {
u1 += pro[i] * i; // C1的平均值
}
v = w0 * w1 * pow(u1 - u0, 2);
if (v >delta) {
delta = v;
thresh = T;
}
T += 1;
}
}
Mat OTSU = Mat::zeros(Size(img.cols, img.rows), img.type());
for (int row = 0; row< img.rows; row++) {
for (int col = 0; col< img.cols; col++) {
if (img.at(row, col) >int(thresh)) {
OTSU.at(row, col) = 255;
}
}
}
imshow("OTSU阈值分割", OTSU);
waitKey(0);
destroyAllWindows();
}
void Sobel(Mat img) { // 基于Prewitt算子的阈值分割
Mat Sobel_City = Mat::zeros(Size(img.cols, img.rows), img.type()); //用于计算城市距离的空白图像
Mat Sobel_Ojld = Mat::zeros(Size(img.cols, img.rows), img.type()); //用于计算欧几里得距离的空白图像
for (int row = 1; row< img.rows - 1; row++) {
for (int col = 1; col< img.cols - 1; col++) { // 由于像素之间的加减可能会小于零，因此记得加上绝对值函数fabs()
Sobel_City.at(row, col) = saturate_cast(fabs(img.at(row - 1, col + 1) - img.at(row - 1, col - 1) + 2*img.at(row, col + 1) - 2*img.at(row, col - 1) + img.at(row + 1, col + 1) - img.at(row + 1, col - 1)) + fabs(img.at(row + 1, col - 1) - img.at(row - 1, col - 1) + 2*img.at(row + 1, col) - 2*img.at(row - 1, col) + img.at(row + 1, col + 1) - img.at(row - 1, col + 1)));
}
}
for (int row = 1; row< img.rows - 1; row++) {
for (int col = 1; col< img.cols - 1; col++) {
Sobel_Ojld.at(row, col) = saturate_cast(sqrt(pow(img.at(row - 1, col + 1) - img.at(row - 1, col - 1) + 2*img.at(row, col + 1) - 2*img.at(row, col - 1) + img.at(row + 1, col + 1) - img.at(row + 1, col - 1), 2) + pow(img.at(row + 1, col - 1) - img.at(row - 1, col - 1) + 2*img.at(row + 1, col) - 2*img.at(row - 1, col) + img.at(row + 1, col + 1) - img.at(row - 1, col + 1), 2)));
}
}
imshow("Sobel图像(街区距离)", Sobel_City);
imshow("Sobel图像(欧几里得距离)", Sobel_Ojld);
int pixel_num01[256] = { 0 }; //数组记得初始化为0，即未统计数量之前各个像素的个数都是0
for (int row = 0; row< img.rows; row++) {
for (int col = 0; col< img.cols; col++) {
pixel_num01[Sobel_Ojld.at(row, col)] += 1; //遍历到的像素值作为索引，次数+1
}
}
int times = 0; //定义遍历数组的变量，从0开始，依次输出0到255每个像素值的数目
while (times<= 255) {
cout<< "像素值"<< times<< "的数目为： "<< pixel_num01[times]<< endl; // 遍历输出
times++; // 不要忘了自增
}
//得到10作为分割阈值
for (int row = 0; row< img.rows; row++) {
for (int col = 0; col< img.cols; col++) { //遍历所有像素点进行判断
if (Sobel_Ojld.at(row, col)< 70) {
Sobel_Ojld.at(row, col) = 0; // 小于阈值赋值为0，否则赋值255
}
else {
Sobel_Ojld.at(row, col) = 255;
}
}
}
imshow("欧几里得阈值分割", Sobel_Ojld);
int pixel_num02[256] = { 0 }; //数组记得初始化为0，即未统计数量之前各个像素的个数都是0
for (int row = 0; row< img.rows; row++) {
for (int col = 0; col< img.cols; col++) {
pixel_num01[Sobel_City.at(row, col)] += 1; //遍历到的像素值作为索引，次数+1
}
}
int time = 0; //定义遍历数组的变量，从0开始，依次输出0到255每个像素值的数目
while (times<= 255) {
cout<< "像素值"<< time<< "的数目为： "<< pixel_num02[time]<< endl; // 遍历输出
time++; // 不要忘了自增
}
OTSU(Sobel_Ojld, pixel_num01); // OTSU阈值分割
OTSU(Sobel_City, pixel_num02);
}
你是否还在寻找稳定的海外服务器提供商？创新互联www.cdcxhl.cn海外机房具备T级流量清洗系统配攻击溯源，准确流量调度确保服务器高可用性，企业级服务器适合批量采购，新人活动首月15元起，快前往官网查看详情吧
文章标题：C++手敲基于梯度图和像素数量数组的OTSU阈值分割-创新互联
本文网址：http://cqcxhl.cn/article/dioeed.html