全局阈值处理 ostu算法

otsu算法

下面是我写的代码，也会参照别人写的

/**
* @brief       ：直方图 统计各个灰度级的数量
* @param [in]  : src_image 源图像
*                height    高
*                width     宽
* @param [out] : histogram 各个灰度级统计结果
* @retval      ：none
*/
void MyAlgorithm::Cal_Histogram(unsigned char *src_image, int height, int width, int *histogram)
{
	unsigned char *psrc_image = src_image;
	int *phistogram = histogram;
	int img_height = height;
	int img_width  = width;
	int i = 0;
	
	for (i = 0; i < (img_height*img_width); i++)
	{
		phistogram[psrc_image[i]]++;
	}
}

/**
* @brief       ：计算各个灰度级的概率
* @param [in]  : histogram 
*                height    高
*                width     宽
* @param [out] : hist_prob 各个灰度级概率
* @retval      ：none
*/
void MyAlgorithm::Cal_Hist_Prob(int *histogram, int height, int width,double *hist_prob)
{
	int *phistogram = histogram;
	double *phist_prob = hist_prob;
	int img_height = height;
	int img_width = width;
	int i = 0;

	for (i = 0; i < 256; i++)
	{
		hist_prob[i] = (double)histogram[i] / (double)(img_height * img_width);
	}
}

/**
* @brief       ：计算各个灰度级的概率累加和
* @param [in]  : hist_prob 灰度级概率
* @param [out] : none
* @retval      ：灰度级概率累加值
*/
double MyAlgorithm::Cal_Hist_Prob_Sum(double *hist_prob, int num)
{
	int i = 0;
	double sum = 0.0f;
	double *phist_prob = hist_prob;

	for (i = 0; i < num; i++)
	{
		sum += phist_prob[i];
	}
	return sum;
}

/**
* @brief       ：计算灰度级累加均值
* @param [in]  : hist_prob 灰度级概率
* @param [out] : hist_prob 各个灰度级概率
* @retval      ：none
*/
double MyAlgorithm::Cal_Ave_Sum(double *hist_prob, int num)
{
	double *phist_prob = hist_prob;
	double sum = 0.0f;
	int i = 0;

	for (i = 0; i < num; i++)
	{
		sum += i*phist_prob[i];
	}
	return sum;
}

/**
* @brief       ：计算类间方差
*                ((mg * p1(k) - m1(k))*(mg * p1(k) - m1(k)))/(p1(k)(1-p1(k)))
* @param [in]  : hist_prob 灰度级概率
* @param [out] : none
* @retval      ：灰度级概率累加值
*/
double MyAlgorithm::Cal_Sigma(double *hist_prob, double mg, int num)
{
	double p1 = 0.0;
	double m1 = 0.0;
	double sigma = 0.0;

	p1 = Cal_Hist_Prob_Sum(hist_prob, num);

	m1 = Cal_Ave_Sum(hist_prob,num);

	if ((fabs(p1) < 1e-15) || (fabs(1 - p1) < 1e-15))
		return 0.0;
	else
	{
		sigma = pow(mg*p1-m1,2)/(p1*(1-p1));
		return sigma;
	}
}

/**
* @brief       ：otsu阈值处理
* @param [in]  : src_image 源图像
*                height    高
*                width     宽
* @param [out] : dst_iamge 目标图像
* @retval      ：none
*/
void MyAlgorithm::Get_Otsu(unsigned char *src_image, unsigned char *dst_iamge, int height, int width)
{
	int i, j;

	unsigned char *psrc_image = src_image;
	unsigned char *pdst_image = dst_iamge;
	int img_height = height;
	int img_width = width;
	int index = 0;

	int histogram[256] = {0};                  //统计各个灰度级数量数组
	double hist_prob[256] = {0.0};            //各个灰度级概率数组
	double hist_prob_sum[256] = {0.0};        //灰度级累加和概率
	double ave_sum[256] = {0.0};              //累加均值

	double sigma[256] = { 0.0 };
	double mg = 0.0;
	double max_sigma = 0.0;
	int max_sigma_count = 0;

	int max_simga_ave = 0;
	int temp = 0;

	Cal_Histogram(psrc_image, img_height, img_width, histogram);
	Cal_Hist_Prob(histogram, img_height, img_width, hist_prob);
	mg = Cal_Ave_Sum(hist_prob,256);

	for (i = 0; i < 256; i++)
		sigma[i] = Cal_Sigma(hist_prob, mg, i);

	max_sigma = sigma[0];

	for (i = 1; i < 256; i++)
	{
		if (fabs(max_sigma - sigma[i]) < 1e-15)
		{
			temp += i;
			max_sigma_count++;
		}
		else if (sigma[i] > max_sigma)
		{
			max_sigma = sigma[i];
			max_sigma_count = 1;
			temp = i;
		}
	}
	max_simga_ave = temp / max_sigma_count;

	for (i = 0; i < height; i++)
	{
		for (j = 0; j < width; j++)
		{
			index = i * width + j;
			if (psrc_image[index] > max_simga_ave)
				pdst_image[index] = 0xff;
			else 
				pdst_image[index] = 0;
		}
	}

}