【数字图像处理】(一)二值化

全局二值法

1.迭代法

迭代法
迭代法是基于逼近的思想,其步骤如下:
(1)求出图象的最大灰度值和最小灰度值,分别记为Rmax和Rmin,令阈值T=(Rmax+Rmin)/2。
(2)根据阈值T将图象的平均灰度值分成两组R1和R2。
(3)分别求出两组的平均灰度值μ1和μ2。
(4)求出新阈值T=(μ1+μ2)/2。

 

2.最大类间方差法OTSU

 

#include "MyIPL.h"

//手动设定阈值
void thresholdSimple(Mat &src,Mat &dst,int threshold )
{
	for (int i=0;i<src.rows;i++)
	{
		for (int j=0;j<src.cols;j++)
		{
			dst.at<uchar>(i,j) = (src.at<uchar>(i,j)>threshold)?0:255;
		}
	}
}


	/*======================================================================*/
	/* 迭代法*/
	/*======================================================================*/
	// nMaxIter：最大迭代次数；nDiffRec：使用给定阀值确定的亮区与暗区平均灰度差异值;
int IterThreshold(Mat img, int nMaxIter, int& iDiffRec)  //阀值分割：迭代法
{
	//图像信息
	int height = img.rows;
	int width = img.cols;
	//int step = width;
	//uchar *data = (uchar*)img->imageData;

	iDiffRec =0;
	int F[256]={ 0 }; //直方图数组
	int iTotalGray=0;//灰度值和
	int iTotalPixel =0;//像素数和
	uchar bt;//某点的像素值

	uchar iThrehold,iNewThrehold;//阀值、新阀值
	uchar iMaxGrayValue=0,iMinGrayValue=255;//原图像中的最大灰度值和最小灰度值
	uchar iMeanGrayValue1,iMeanGrayValue2;

	//获取(i,j)的值，存于直方图数组F
	for(int i=0;i<width;i++)
	{
		for(int j=0;j<height;j++)
		{
			bt = img.at<uchar>(j,i);
			if(bt<iMinGrayValue)
				iMinGrayValue = bt;
			if(bt>iMaxGrayValue)
				iMaxGrayValue = bt;
			F[bt]++;
		}
	}

	iThrehold =0;//
	iNewThrehold = (iMinGrayValue+iMaxGrayValue)/2;//初始阀值
	iDiffRec = iMaxGrayValue - iMinGrayValue;

	for(int a=0;(abs(iThrehold-iNewThrehold)>0.5)&&a<nMaxIter;a++)//迭代中止条件
	{
		iThrehold = iNewThrehold;
		//小于当前阀值部分的平均灰度值
		for(int i=iMinGrayValue;i<iThrehold;i++)
		{
			iTotalGray += F[i]*i;//F[]存储图像信息
			iTotalPixel += F[i];
		}
		iMeanGrayValue1 = (uchar)(iTotalGray/iTotalPixel);
		//大于当前阀值部分的平均灰度值
		iTotalPixel =0;
		iTotalGray =0;
		for(int j=iThrehold+1;j<iMaxGrayValue;j++)
		{
			iTotalGray += F[j]*j;//F[]存储图像信息
			iTotalPixel += F[j];    
		}
		iMeanGrayValue2 = (uchar)(iTotalGray/iTotalPixel);

		iNewThrehold = (iMeanGrayValue2+iMeanGrayValue1)/2;        //新阀值
		iDiffRec = abs(iMeanGrayValue2 - iMeanGrayValue1);
	}

	//cout<<"The Threshold of this Image in imgIteration is:"<<iThrehold<<endl;
	return iThrehold;
}


	/*======================================================================*/
	/* OTSU global thresholding routine */
	/* takes a 2D unsigned char array pointer, number of rows, and */
	/* number of cols in the array. returns the value of the threshold */
	/*parameter: 
	*image --- buffer for image
	rows, cols --- size of image
	x0, y0, dx, dy --- region of vector used for computing threshold
	vvv --- debug option, is 0, no debug information outputed
	*/
	/*
	OTSU 算法可以说是自适应计算单阈值（用来转换灰度图像为二值图像）的简单高效方法。
	下面的代码最早由 Ryan Dibble提供，此后经过多人Joerg.Schulenburg, R.Z.Liu 等修改，补正。
	算法对输入的灰度图像的直方图进行分析，将直方图分成两个部分，使得两部分之间的距离最大。
	划分点就是求得的阈值。
	*/
	/*======================================================================*/
int otsu(Mat image)
{
	int w = image.cols;
	int h = image.rows;

	unsigned char*np; // 图像指针
	unsigned char pixel;
	int thresholdValue=1; // 阈值
	int ihist[256]; // 图像直方图，256个点

	int i, j, k; // various counters
	int n, n1, n2, gmin, gmax;
	double m1, m2, sum, csum, fmax, sb;

	// 对直方图置零...
	memset(ihist, 0, sizeof(ihist));

	gmin=255; gmax=0;
	// 生成直方图
	for (i =0; i < h; i++) 
	{
		for (j =0; j < w; j++) 
		{
			pixel = image.at<uchar>(i,j);
			ihist[ pixel]++;
			if(pixel > gmax) gmax= pixel;
			if(pixel < gmin) gmin= pixel;
		}
	}

	// set up everything
	sum = csum =0.0;
	n =0;

	for (k =0; k <=255; k++) 
	{
		sum += k * ihist[k]; /* x*f(x) 质量矩*/
		n += ihist[k]; /* f(x) 质量 */
	}

	if (!n) 
	{
		// if n has no value, there is problems...
		//fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");
		thresholdValue =160;
		goto L;
	}

	// do the otsu global thresholding method
	fmax =-1.0;
	n1 =0;
	for (k =0; k <255; k++) 
	{
		n1 += ihist[k];
		if (!n1) { continue; }
		n2 = n - n1;
		if (n2 ==0) { break; }
		csum += k *ihist[k];
		m1 = csum / n1;
		m2 = (sum - csum) / n2;
		sb = n1 * n2 *(m1 - m2) * (m1 - m2);
		/* bbg: note: can be optimized. */
		if (sb > fmax)
		{
			fmax = sb;
			thresholdValue = k;
		}
	}

L:
	for (i =0; i < h; i++) 
	{
		//np = (unsigned char*)(image->imageData + image->widthStep*i);
		for (j =0; j < w; j++) 
		{
			if(image.at<uchar>(i,j) >= thresholdValue)
				np[j] =255;
			else np[j] =0;
		}
	}

	//cout<<"The Threshold of this Image in Otsu is:"<<thresholdValue<<endl;
	return(thresholdValue);
}