【图像处理】灰度加权距离变换（GWDT）

 

【fishing-pan：https://blog.csdn.net/u013921430转载请注明出处】

前言

       距离变换（distance transform，DT）在图像处理、计算机视觉等领域有非常多的运用，例如骨架化、目标细化、黏连目标分离等。在Matlab和OpenCV中都有进行距离变换的函数。Matlab中有bwdist函数用于二值化图像的距离变换， OpenCV中有cvDistTransform函数用于进行二值化图像的距离变换。

      除此之外，Matlab中还提供了graydist函数用于灰度加权距离变换，但是graydist居然要输入一个种子点，这让我觉得很神奇，因为在我的印象中灰度加权距离变换是不需要种子点的！自己测试运行后就发现，这是扯淡呢！算出来的结果是图像中各个点到种子点的距离加权变换。所以，我必须来讲讲我所理解的不需要种子点的灰度加权距离变换。

灰度加权距离变换

       一般的距离变换需要首先对图像进行二值化，这样的过程会损失掉图像中某些有用的信息，因为单纯的距离变换获得的是前景点到背景的最短距离，损失了灰度信息所代表的物理意义。因此在距离变换的基础上加上图像灰度信息，就有了灰度加权距离变换——GWDT（gray-weighted image distance transform），也可以称之为灰度尺度距离变换——GSDT（gray-scale image distance transform）。

       理解了距离变换，灰度加权距离变换就很好理解了，就是在距离变换的基础上加上了灰度信息。那么如何操作呢？我采用类似于区域增长的方法，在快速行进法的框架下来进行实现它，基本步骤如下；

1.    设定阈值分离前景和背景，将背景点的状态设置为ALIVE（其实所有的背景点都是seed point），将前景设为FAR，前景点的值初始化为无限大；

2.    初始化边缘点，在FAR点中选取8邻域（对二维图像进行处理）中有背景点的点，这些点就是边缘点，将边缘点的值赋予原图中的值，更改其状态为TRAIL，并将其放入一个队列Q中；

3.    在队列Q中不放回地取出点，假设取出的点记为y，计算y周围的8个点x，根据以下公式计算x，如果x的值被更新了，判断其是否已经在队列中，如果不在，那么就将其放入队列中。（在我实现的时候，考虑到效率的因素，并没有判断其是否已经在队列中，而是将其直接放入队列中，这个跟queue类和C++函数传参有关系）

4.    循环第三步，直到队列Q为空。

       从上面的步骤中可以看出，实现过程与快速行进算法还是存在一定的差异的，如果按照快速行进算法的思想，那么对于点y，只会计算周围8临域中属性是FAR的点的值，不会对ALIVE的点进行更新。但是上述的过程，会对ALIVE点进行更新。我这样做是为了尽量保证每个点在变换过后值尽量小，大家可以自行修改代码来决定要不要对ALIVE点进行更新。 

---

       如下图所示，假设图中蓝色的点为初始化后的边缘点，他的状态为TRAIL，点A的状态为FAR；依次遍历S1-S4周围的8个点，很显然A是这四个点的共同邻域，A取四个点对其进行变换后的最小值；

 

                                                                         图 1 迭代过程示意图标题标题

 

处理结果

                                                                                        图 2 原始图像

 

 

                                                                                  图 3  处理后的图像

---

结果分析

       原图中我加入了一些噪声，在GWDT后，前景点特别是远离背景的前景点的灰度值相对被增强了，这样背景点就变相的被抑制了。此外，通过设定阈值可以分离区域内的黏连的大米，从而自动进行记数。等等等等，GWDT的应用还有很多。总的来说，GWDT是在传统的距离变换的基础上保留和利用了灰度信息，更加具有物理意义和价值。

代码

       为了方便大家学习，为大家提供源码；

//-------------------------
//潘正宇 ，2018.05.10
//灰度尺度距离变换，GSDT/GWDT
//-------------------------

#include <iostream>
#include <vector>
#include <opencv2/opencv.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <queue>
#include <string>

using namespace std;
using namespace cv;


#define INF 3.4e+38                        // 无穷大的点
enum{ ALIVE = -1, TRIAL = 0, FAR = 1 };    //设定三种状态；


bool find(queue<int> que,int value)
{
	int size = que.size();
	for (int i = 0; i <size; ++i)
	{
		int a = que.front();
		if (a == value)
		{
			return true;
		}
		que.pop();
	}
	return false;
}


bool GSDT(string &picPath)
{
	Mat Scr = imread(picPath, 0);    //以灰度图的模式读入图像

	imshow("原始图", Scr);

	Mat gsdtImg;                     //创建GSDT后的图
	gsdtImg.create(Scr.size(), CV_32FC1);

	int row = Scr.rows;
	int col = Scr.cols;

	Mat  mat_mean, mat_stddev;	
	meanStdDev(Scr, mat_mean, mat_stddev);
	double mean;
	mean = mat_mean.at<double>(0, 0);   //获取图像均值

	int *state= new int[col*row];

	for (int i = 0; i < row; i++)
	{
		for (int j = 0; j < col; j++)
		{
			if (Scr.at<uchar>(i, j)<(mean))
			{
				gsdtImg.at<float>(i, j) = Scr.at<uchar>(i, j);       //背景点是ALIVE的
				state[i*col + j] = ALIVE;
			}
			else
			{
				gsdtImg.at<float>(i, j) = INF;
				state[i*col + j] = FAR;
			}

		}
	}
	
	queue<int> TrailQue;                     //定义队列用于存放TRAIL的点
	for (int i = 1; i < row-1; i++)
	{
		for (int j = 1; j < col-1; j++)
		{
			if (state[i*col + j] == FAR)     //如果这个点是前景点，那么搜寻这个点周围是否有背景点，如果是那么它就是边缘点；
			{
				for (int o = -1; o <= 1; o++)
				{
					for (int p = -1; p <= 1; p++)
					{
						if (state[(i + o)*col + j + p] == ALIVE)        //找到所有的边缘点，并且将其放入队列TrailQue；
						{
							state[i*col + j] = TRIAL;                   
							gsdtImg.at<float>(i, j) = Scr.at<uchar>(i, j);
							TrailQue.push(i*col + j);
							break;
						}
					}
				}
			}
			
		}
	}
	

	while (!TrailQue.empty())
	{
		int P_row = TrailQue.front() / col;    ///获取TrailQue中点的坐标
		int P_col = TrailQue.front() % col;

		if (P_row < 1){ P_row = 1; }if (P_row>row - 2){ P_row = row - 2; }
		if (P_col < 1){ P_col = 1; }if (P_col>col - 2){ P_col = col - 2; }

		for (int o = -1; o <= 1; o++)
		{
			for (int p = -1; p <= 1; p++)
			{
				int N = (P_row + o)*col + P_col + p;

				double len = sqrt(o*o + p*p);
				float gs = gsdtImg.at<float>(P_row, P_col) + Scr.at<uchar>(P_row + o, P_col + p)*len;

				//---比较该点现有的GSDT值与P点给它的值；
				if (gsdtImg.at<float>(P_row + o, P_col + p) > gs)
				{
					state[N] = TRIAL;
					gsdtImg.at<float>(P_row + o, P_col + p) = gs;

					TrailQue.push(N);

					/*if (!find(TrailQue, N))
					{
						TrailQue.push(N);
					}*/
					
				}
					
				
			}
		}
		

		TrailQue.pop();                        //从TrailQue中删除这个点
	}

	//寻找最大值；
	float Max = 0;
	for (int i = 0; i < row; i++)
	{
		for (int j = 0; j < col; j++)
		{
			if (gsdtImg.at<float>(i, j)>(INF / 100))
			{
				gsdtImg.at<float>(i, j) = 0;
			}
			if (gsdtImg.at<float>(i, j)>Max)//&&gsdtImg.at<float>(i, j)<(INF/10)
			{
				Max = gsdtImg.at<float>(i, j);
			}
			
		}
	}

	//将图像像素区间压缩到0-255
	Mat Dst;
	Dst = Scr.clone();
	for (int i = 0; i < row; i++)
	{
		for (int j = 0; j < col; j++)
		{
			Dst.at<uchar>(i, j) = 255 * gsdtImg.at<float>(i, j) / (Max+1);
		}
	}

	imwrite("temp.bmp", Dst);
	
	//归一化
	gsdtImg = gsdtImg / (Max+1);
	
	
	cout << gsdtImg.at<float>(20, 206);
	imshow("GSDT",gsdtImg);
	waitKey(0);

	return true;
}


void main()
{
	string picPath = "G:\\博客\\图像处理\\灰度尺度距离变换\\rice.jpg";
	GSDT(picPath);

	system("pause");
	return;
}