Exercise5--阈值分割

由于阈值处理直观，实现简单且计算速度快，因此阈值处理再图像分割应用处于核心地位

---

假设一个图像的灰度直方图对应于图像的$f(x,y)$,而该图像由暗色背景上的较亮物体组成，那么，就可以设置一个全局阈值T，然后对于每一个$f(x,y)>T$的点称之为对象点，其余的称之为背景点，从而达到了分割图像的目的。

$$g(x,y)=\begin{cases}1,&f(x,y)>T\\0,&f(x,y)\leq T\end{cases}$$

基本的全局阈值处理

当物体的背景像素的灰度分布十分明显时，可以用适用于整个图像的单个（全局）阈值。通常图像之间有较大的变化，使用手动设置的全局阈值是一种方法，当然也有着对图像进行自动阈值估算的算法

迭代法

1，为全局阈值$T_{0}$选择一个初始估计值(一般为中位数)

2，使用 $T_{0}$把图像分割成两个部分 $R_{1}$和$R_{2}$,并计算其平均灰度 $m_{1}$和$m_{2}$

$$m_{1} = \frac{\sum_{i=0}^{T_{0}}i*n_{i}}{\sum_{i=0}^{T_{0}}n_{i}},m_{2} = \frac{\sum_{i=T_{0}}^{L-1}i*n_{i}}{\sum_{i=T_{0}}^{L-1}n_{i}}$$

3,计算新阈值

$$T_{0}=\frac{m_{1}+m_{2}}{2}$$

4,重复步骤2到4，直到$T_{0} <\Delta T$

这样就可以计算出一个比较合适的阈值$T_{0}$

实现

void thresholdByIterate(const ImageUtil::ImageData& data)
{
	ImageUtil::GrayHistogram histogram = ImageUtil::getHistogram(data);
	histogram.normalize();
    //开始阈值设置为中位数，因此需要计算
	std::vector<int> his(data.width * data.height);
	for (int i = 0; i < data.width * data.height; i++)
	{
		his[i] = data.pImg[i];
	}

	std::sort(his.begin(), his.end());
	double t0 = 0,t1 = his[his.size() / 2];


    //不断计算阈值左右的平均灰度，从而达到左右比较平衡的目的
	while (std::abs(t0 - t1) > 10)
	{
		double a0 = 0,n0 = 0, a1 = 0,n1 = 0;
		for (int i = 0; i < t1; i++)
		{
			a0 += histogram.gray[i] * i;
		}
		for (int i = 0; i < t1; i++)
		{
			n0 += histogram.gray[i];
		}

		for (int i = t1; i < 255; i++)
		{
			a1 += histogram.gray[i] * i;
		}
		for (int i = t1; i < 255; i++)
		{
			n1 += histogram.gray[i];
		}

		t0 = t1;
		t1 = (a0 / n0 + a1 / n1) * 0.5;
	}

	//根据阈值重新设置图像
	ImageUtil::ImageData img = data;
	BYTE *imgData = new BYTE[data.width * data.height];
	for (int i = 0; i < data.width * data.height; i++)
	{
		imgData[i] = data.pImg[i] > t1 ? 1 : 0;
	}

    //输出黑白图
	std::string name("bitmap/threshold_by_iterate");
	name.append("_")
		.append(std::to_string(t1))
		.append(".bmp");


	img.pImg = imgData;
	ImageUtil::outputBlackWhiteImage(img, name);

	delete[] imgData;
}

Otsu法

阈值处理可以视为一种统计决策理论问题。其目的在于把像素分为两个或多个组，使得分组的时候引入的平均误差最小。

而Otsu法则是一种很好的方案，他的思想在于，使得分组之间的类间方差最大化，这么，也就是相当于每一个类里面的方差是最小的，从而达成了一个比较好的分割的目的。

Otsu方法还有一个特性，就是他可以完全在一个图片上的直方图进行运算。

由于Otsu方法的思想在于取得最大的类间方差，那么，我们需要做的就是计算出每一个类里的平均灰度值$m$，以及他们之间的方差即可。

证明

以下证明来自于冈萨雷斯的《数字图像处理（第三版）》

首先，我们选择一个阈值$T(k)=k,(0<k<L-1)$,然后把输入的图像分为两类$C_{1}$和$C_{2}$,那么$C_{1}$当中所有的所有的点由所有灰度值在$[0,k]$内的像素组成，而$C_{2}$同理。

那么我们可以得知，类$C_{1}$的概率$P_{1}$可以这样得到

$$P_{1}(k)=\sum_{i=0}^kP_{i}$$

类$C_{2}$同理

$$P_{2}(k)=\sum_{i=k+1}^{L-1}P_{i}=1-P_{1}(k)$$

而分配到类$C_{1}$的平均灰度值为

$$m_1(k)=\sum_{i=0}^kiP(i/C_1)=\frac{1}{P_1(k)}\sum_{i=0}^kiP_i$$

同样$C_{2}$同理。

而我们同样可以根据直方图计算出全局的平均灰度值

$$m_G=\sum_{i=0}^{L-1}iP_i$$

我们可以这样验证这个公式

$$\begin{cases}P_1m_1+P_2m_2&=m_G\\P_1+P_2&=1\end{cases}$$

在拥有了每一个类的平均灰度值与全局灰度值之后，我们就可以计算类间的方差了。

$$\begin{align}\sigma_{b}^2&=P_1(m_1-m_G)^2+P_2(m_2-m_G)^2 \\&=P_1P_2(m_1-m_2)^2\\&=\frac{(m_GP_1-m)^2}{P_1(1-P_1)}\end{align}$$

由于$m_G$只需要计算一次，所以，我们只需要对应每一个选定的阈值k都计算其P_1，m_k既可以计算出所有阈值下的方差了。

最后，我们只需要寻找到方差最大的一个阈值，设置为最终阈值即可。(若有多个相同大小的最大阈值，就取其平均值)。

实现

根据公式，即可以写出其中的算式实现。

//计算方差
double otsuVariance(const int k, const double mG,const ImageUtil::GrayHistogram& histogram)
{
	return std::pow(mG*otsuP(k, histogram) - otsuM(k, histogram), 2) / (otsuP(k, histogram) * (1 - otsuP(k, histogram)));
}

//计算平均灰度值
double otsuM(const int k,const ImageUtil::GrayHistogram& histogram)
{
	double result = 0;
	for(int i = 0;i < k;i++)
	{
		result += i * histogram.gray[i];
	}

	return result;
}

//计算概率P
double otsuP(const int k,const ImageUtil::GrayHistogram& histogram)
{
	double result = 0;
	for(int i = 0;i < k;i++)
	{
		result += histogram.gray[i];
	}

	return result;
}

然后，我们就可以开始利用这些公式计算了。

1，输入归一化的直方图

2，计算出全局灰度

3，对于每一个灰度分量阈值k都计算出其类间方差

4，选出其中的方差最大的作为阈值k，若有相同的，取平均值

从而既可以实现Ostu方法了。

void otsu(const ImageUtil::ImageData& data)
{
	ImageUtil::GrayHistogram histogram = ImageUtil::getHistogram(data);
	histogram.normalize();

	const int len = 255;

    //全局阈值
	const double mG = otsuM(len, histogram);
	double delta[len];
    //计算每一个灰度分量的阈值k
	for (int i = 0; i < len; i++)
	{
		delta[i] = otsuVariance(i, mG, histogram);
	}

	double max = -1;
	for (double i : delta)
	{
		if(i > max)
		{
			max = i;
		}
	}

	std::vector<int> maxList;
	for(int i = 0;i < len;i++)
	{
		if (delta[i] == max)
			maxList.push_back(i);
	}

	int k = 0;
	for (int i : maxList)
	{
		k += i;
	}

	k /= maxList.size();

    //根据阈值重新设置图像
	ImageUtil::ImageData img = data;
	BYTE *imgData = new BYTE[data.width * data.height];
	for (int i = 0; i < data.width * data.height; i++)
	{
		imgData[i] = data.pImg[i] > k ? 1 : 0;
	}

    //输出黑白图
	std::string name("bitmap/threshold_by_otsu");
	name.append("_")
		.append(std::to_string(k))
		.append(".bmp");

	img.pImg = imgData;
	ImageUtil::outputBlackWhiteImage(img, name);
}

基于边缘的灰度分割

由于直方图很少依赖物体与背景的大小，因此，有时候会出现阈值处理失败的情况。比如，一个大的背景当中有很多噪点，这样会污染了直方图。

因此，在直方图的处理上有一种改进就是利用图像的边缘。因为图像的边缘四周是背景或者是前景图的概率是很高的。

这种算法的思想在于

1，计算图像的边缘(拉普拉斯算子之类的)

2，设置阈值T，将边缘图当中所有大于T的像素点设为1，否则为0

3，根据像素图当中1的点的位置，只使用这些位置的点，计算出原图的直方图

4，对直方图进行全局分割（例如Otsu方法）

Laplace + Otsu 实现

void laplaceOstu(const ImageUtil::IMGDATA& data)
{
	BYTE * newData = new BYTE[data.length];

    //利用Laplace算法计算出边缘
	for (int i = 0; i < data.height; i++)
	{
		for (int j = 0; j < data.width; j++)
		{
			int up, down, left, right;
			if (i == 0)
				up = 0;
			else
				up = i - 1;

			if (i == data.height - 1)
				down = data.height - 1;
			else
				down = i + 1;

			if (j == 0)
				left = 0;
			else
				left = j - 1;

			if (j == data.width - 1)
				right = data.width - 1;
			else
				right = j + 1;

			newData[i * data.width + j] = ImageUtil::clamp(
				1 * data.pImg[up * data.width + left] + 1 * data.pImg[up * data.width + j] + 1 * data.pImg[up * data.width + right] +
				1 * data.pImg[i * data.width + left] + -8 * data.pImg[i * data.width + j] + 1 * data.pImg[i * data.width + right] +
				1 * data.pImg[down * data.width + left] + 1 * data.pImg[down * data.width + j] + 1 * data.pImg[down * data.width + right]);
		}

	}

    //这里设置的阈值T为150
	for(int i = 0;i < data.width * data.height;i++)
	{
		newData[i] = newData[i] > 150 ? 1 : 0;
	}

	ImageUtil::IMGDATA imgData = data;
	imgData.pImg = newData;

    //计算出灰度图
	ImageUtil::GRAYHISTOGRAM grayhistogram;
	grayhistogram.pixelCount = 0;
	int point = -1;
	for (int i = 0; i < data.height; i++)
	{
		for (int j = 0; j < data.width; j++)
		{
			
			if(newData[++point] == 1)
			{
				grayhistogram.gray[data.pImg[point]]++;
				grayhistogram.pixelCount++;
			}
			
		}
	}
    
    //输出灰度图
	grayhistogram.normalize();
	ImageUtil::outputHistogram(grayhistogram, "bitmap/laplace_histogram.bmp");
	const int len = 255;

    //利用Ostu方法计算出阈值
	const double mG = otsuM(len, grayhistogram);
	double delta[len];
	for (int i = 0; i < len; i++)
	{
		delta[i] = otsuVariance(i, mG, grayhistogram);
	}

	double max = -1;
	for (double i : delta)
	{
		if (i > max)
		{
			max = i;
		}
	}

	std::vector<int> maxList;
	for (int i = 0; i < len; i++)
	{
		if (delta[i] == max)
			maxList.push_back(i);
	}

	int k = 0;
	for (int i : maxList)
	{
		k += i;
	}

	k /= maxList.size();

	ImageUtil::ImageData img = data;
	BYTE *mg = new BYTE[data.width * data.height];
	for (int i = 0; i < data.width * data.height; i++)
	{
		mg[i] = data.pImg[i] > k ? 1 : 0;
	}
	
    //输出最终的黑白图
	std::string name("bitmap/threshold_by_laplace_otsu");
	name.append("_")
		.append(std::to_string(k))
		.append(".bmp");

	img.pImg = mg;
	ImageUtil::outputBlackWhiteImage(img, name);

	delete[] newData;
	delete[] mg;

}

完整源码