C++ opencv基于OTSU图像多阈值分割

vector<double> sub(vector<int> A, vector<double> B) //减法
{
	vector<double> C;
	C.resize(3);
	for (int i = 0; i < A.size(); i++)
	{
		C[i] = A[i] - B[i];
	}
	return C;
}

vector<double> multi(vector<double> A, int b)//乘法
{
	vector<double> C;
	C.resize(3);
	for (int i = 0; i < A.size(); i++)   
	{
		C[i] = A[i] * b;
	}
	return C;
}

vector<double> multi(vector<double> A, vector<double> B)//乘法
{
	vector<double> C;
	C.resize(3);
	for (int i = 0; i < A.size(); i++)
	{
		C[i] = A[i] * B[i];
	}
	return C;
}

vector<double> divi(vector<double> A, vector<double> B)//处罚
{
	vector<double> C;
	C.resize(3);
	for (int i = 0; i < A.size(); i++)
	{
		if(B[i]==0)
		{
			continue;
		}
		C[i] = A[i] / B[i];
	}
	return C;
}



/*
* 实现python的linspace函数功能
*/
std::vector<double> linspace(int start_in, int end_in, int num_in)
{

	std::vector<double> linspaced;

	double start = static_cast<double>(start_in);
	double end = static_cast<double>(end_in);
	double num = static_cast<double>(num_in);

	if (num == 0) { return linspaced; }
	if (num == 1)
	{
		linspaced.push_back(start);
		return linspaced;
	}

	double delta = (end - start) / (num - 1);

	for (int i = 0; i < num - 1; ++i)
	{
		linspaced.push_back(start + delta * i);
	}
	linspaced.push_back(end); 
							
	return linspaced;
}
/*
* vector容器第n个到m个元素求和
*/
double compute_sum(vector<double>& Vec, int start, int end)
{
	double total = 0;
	for (int i = start; i < end; i++)
	{
		total += Vec[i];
	}	
	return total;
}

/*
* vector容器dot点乘
*/
double compute_dot(vector<double>& Vec1, vector<double>& Vec2, int start, int end)
{
	double total = 0;
	double temp = 0;
	for (int i = start; i < end; i++)
	{
		temp = Vec1[i]*Vec2[i];
		total = total + temp;
	}
	return total;
}

void otsu3layers(Mat &input)
{
	//Mat img = input.clone();
	int high = input.rows;
	int wide = input.cols;
	Mat hist;
	const int bins[1] = { 256 };
	float hranges[2] = { 0,256 };
	const float* ranges[1] = { hranges };
	calcHist(&input, 1, 0, Mat(), hist, 1, bins, ranges);//hist256*1，每个灰度级中所含的像素数
	vector<double>p;//每个像素所占总像素数的比例
	std::vector <double>lsp;
	double temp = 0;
	for (int i = 0; i < hist.rows; i++)
	{
		for (int j = 0; j < hist.cols; j++)
		{
			temp = hist.at<float>(i, j) / (high * wide);
			p.push_back(temp);
		}
	}
	lsp = linspace(1, 256, 256);//.reshape没写
	double maxvar = 0;
	int th1 = 0;
	int th2 = 0;
	int cnt = 0;
	double w0 = 0;
	double w1 = 0;
	double w2 = 0;
	double u0 = 0;
	double u1 = 0;
	double u2 = 0;
	double varbetween = 0;
	vector<double>W;
	vector<int>U;
	vector<double>tempp;

	for (int k1 = 1; k1 < 257; k1++)
	{
		for (int k2 = k1 + 1; k2 < 257; k2++)
		{
			w0 = compute_sum(p, 0, k1);//0.00118
			w1 = compute_sum(p, k1, k2);//0.00076
			w2 = compute_sum(p, k2, 256);//0.998
			u0 = compute_dot(lsp, p, 0, k1);//0.00118
			u1 = compute_dot(lsp, p, k1, k2);//0.00152
			u2 = compute_dot(lsp, p, k2, 256);//53.78
			W.push_back(w0);
			W.push_back(w1);
			W.push_back(w2);
			U.push_back((int)u0);
			U.push_back((int)u1);
			U.push_back((int)u2);
			int UT = accumulate(U.begin(), U.end(), 0);//53
			tempp = multi(W, UT);
			tempp = sub(U, tempp);
			tempp = multi(tempp, tempp);
			tempp = divi(tempp, W);
			varbetween = compute_sum(tempp, 0, 3);
			if (varbetween > maxvar)
			{
				cnt = cnt + 1;
				maxvar = varbetween;
				th1 = k1;
				th2 = k2;
			}
			W.clear();
			U.clear();
			tempp.clear();
			/*cout << "阈值为" << th1 << "和" << th2 << endl;*/
		}
	}
	for (int row = 0; row < high; row++)
	{
		uchar* cur = input.ptr<uchar>(row);
		for (int col = 0; col < wide; col++)
		{
			uchar gray = cur[col];
			/*if (gray < th1)
			{
				gray = 0;
				cur[col] = gray;
			}
			else if (gray > th1 && gray < th2)
			{
				gray = 127;
				cur[col] = gray;
			}
			else
			{
				gray = 255;
				cur[col] = gray;
			}*/
			if (gray < th2)
			{
				gray = 255;
				cur[col] = gray;
			}
			else
			{
				gray = 0;
				cur[col] = gray;
			}
		}
	}

}

参考自：https://www.cnblogs.com/Linkdom/p/16122676.html#%E5%9F%BA%E4%BA%8Eotsu%E7%9A%84%E5%9B%BE%E5%83%8F%E5%A4%9A%E9%98%88%E5%80%BC%E5%88%86%E5%89%B2%E7%AE%97%E6%B3%95