C++手敲基于梯度图和像素数量数组的OTSU阈值分割-创新互联

 一、OTSU算法原理

➢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://cdkjz.cn/article/dioeed.html