itk中的Sobel算子

这篇算是《itk中的特征提取算法（一）》的补丁，Sobel 是一阶微分算子，属于一种常用的图像锐化方法。

作者： Irwin Sobel，1968年提出，两个方向算子：检测水平边缘，检测垂直边缘。

优点：快，简单有效
缺点：不基于灰度对图像进行处理，噪声和条纹都得到了增强。

itk中提供了itkSobelEdgeDetectionImageFilter 和 itkSobelOperator。具体使用方法在参考文献1，2中，而且在《itk中的特征提取算法（一）》中有关于itkSobelEdgeDetectionImageFilter 的代码解释，这里就不再copy了。

下面提供一个自定义的sobel函数：

ImageType *GradientOperator(ImageType *input_data)
{
	ImageType ::PointType origin_temp = input_data->GetOrigin();
	ImageType ::SpacingType spacing_temp = input_data->GetSpacing();
	ImageType ::RegionType inputRegion = input_data->GetLargestPossibleRegion();
	ImageType ::SizeType size = inputRegion.GetSize();
	ImageType ::IndexType start  = inputRegion.GetIndex();
	itk::ImageRegion<2> region(start, size);

	ImageType ::Pointer output = ImageType ::New();
	output_data->SetOrigin(origin_temp);
	output_data->SetSpacing(spacing_temp);
	output_data->SetRegions(region);
	output_data->Allocate();
	output_data->FillBuffer(itk::NumericTraits<PixelType>::Zero);
//
	static const int sizeOfSobelMask = 9;
//Sobel
	static int sobelMaskHor[sizeOfSobelMask] = 
	{
		-1, -2, -1, 
		 0,  0,  0,
		 1,  2,  1
	};

	static int SobelMaskVer[sizeOfSobelMask] = 
	{
		1, 0, -1,
		2, 0, -2,
		1, 0, 1
	};
//Robert~ 所以，只需要把这里的算子换一换就OK啦~不要太简单。
	//static int sobelMaskHor[sizeOfSobelMask] = 
	//{
	//	0,  0, 0, 
	//	0, -1, 0,
	//	0,  0, 1
	//};

	//static int SobelMaskVer[sizeOfSobelMask] = 
	//{
	//	0,  0, 0, 
	//	0,  0, -1,
	//	0,  1, 0
	//};
	//

	typedef itk::ImageRegionIteratorWithIndex<ImageType >   FieldIterator;
	FieldIterator fieldIter( input_data, input_data->GetLargestPossibleRegion());//创建一个像素迭代器
	fieldIter.GoToBegin();//迭代器指向图像数组开始
	ImageType ::IndexType index;

	FieldIterator fieldIter_out( output_data, output_data->GetLargestPossibleRegion());/输出数据像素迭代器
	fieldIter_out.GoToBegin();
	ImageType ::IndexType index_out;

	while( !fieldIter.IsAtEnd()  )//不到最后一个点，一直循环下去
	{
		index=fieldIter.GetIndex();//获得该点坐标
		unsigned char p = 0;//input_data->GetPixel(index);//获得该点像素值

		unsigned char ts = 0, tv = 0,dst = 0;
		ImageType2D::IndexType index_temp;
		
		if((index[0]-1) < 0)
			index[0]++;
		if((index[1]-1) < 0)
			index[1]++;
		if((index[0]+1) > 255)
			index[0]--;
		if((index[1]+1) > 255)
			index[1]--;
//9宫格 = 自己 + 8邻域
		//0
		index_temp[0] = index[0]-1;
		index_temp[1] = index[1]+1;
		p = input_data->GetPixel(index_temp);
		ts += (sobelMaskHor[0] * p );
		tv += (SobelMaskVer[0] * p );
		//1
		index_temp[0] = index[0];
		index_temp[1] = index[1]+1;
		p = input_data->GetPixel(index_temp);
		ts += (sobelMaskHor[1] * p );
		tv += (SobelMaskVer[1] * p );
		//2
		index_temp[0] = index[0]+1;
		index_temp[1] = index[1]+1;
		p = input_data->GetPixel(index_temp);
		ts += (sobelMaskHor[2] * p );
		tv += (SobelMaskVer[2] * p );
		//3
		index_temp[0] = index[0]-1;
		index_temp[1] = index[1];
		p = input_data->GetPixel(index_temp);
		ts += (sobelMaskHor[3] * p );
		tv += (SobelMaskVer[3] * p );
		//4
		index_temp[0] = index[0];
		index_temp[1] = index[1];
		p = input_data->GetPixel(index_temp);
		ts += (sobelMaskHor[4] * p );
		tv += (SobelMaskVer[4] * p );
		//5
		index_temp[0] = index[0]+1;
		index_temp[1] = index[1];
		p = input_data->GetPixel(index_temp);
		ts += (sobelMaskHor[5] * p );
		tv += (SobelMaskVer[5] * p );
		//6
		index_temp[0] = index[0]-1;
		index_temp[1] = index[1]-1;
		p = input_data->GetPixel(index_temp);
		ts += (sobelMaskHor[6] * p );
		tv += (SobelMaskVer[6] * p );
		//7
		index_temp[0] = index[0];
		index_temp[1] = index[1]-1;
		p = input_data->GetPixel(index_temp);
		ts += (sobelMaskHor[7] * p );
		tv += (SobelMaskVer[7] * p );
		//8
		index_temp[0] = index[0]+1;
		index_temp[1] = index[1]-1;
		p = input_data->GetPixel(index_temp);
		ts += (sobelMaskHor[8] * p );
		tv += (SobelMaskVer[8] * p );

		dst = (unsigned char)sqrt( (float)(ts * ts + tv * tv) );
			
		index_out=fieldIter_out.GetIndex();
		output_data->SetPixel(index_out,dst);
		/*
		if(dst > (unsigned char)(0.6180339887 * 255.0) ) //Σ( ° △ °|||)︴ 黄金分割，实际上这里应该用OTUS大致算一个合理的阈值参数
		{
			output_data->SetPixel(index_out,1);
		}
		else 
		{
			output_data->SetPixel(index_out,0);
		}
		*/
		++fieldIter;
		++fieldIter_out;
	} 
	return output_data;
}

分享一个网友写的不错的科普（来自于参考文献3，侵删）：

“首先，我们来开一下计算机是如何检测边缘的。以灰度图像为例，它的理论基础是这样的，如果出现一个边缘，那么图像的灰度就会有一定的变化，为了方便假设由黑渐变为白代表一个边界，那么对其灰度分析，在边缘的灰度函数就是一个一次函数y=kx，对其求一阶导数就是其斜率k，就是说边缘的一阶导数是一个常数，而由于非边缘的一阶导数为零，这样通过求一阶导数就能初步判断图像的边缘了。通常是X方向和Y方向的导数，也就是梯度。理论上计算机就是通过这种方式来获得图像的边缘。

但是，具体应用到图像中你会发现这个导数是求不了的，因为没一个准确的函数让你去求导，而且计算机在求解析解要比求数值解麻烦得多， 所以就想到了一种替代的方式来求导数。就是用一个3×3的窗口来对图像进行近似求导。拿对X方向求导为例，某一点的导数为第三列的元素之和减去第一列元素之和，这样就求得了某一点的近似导数。其实也很好理解为什么它就近似代表导数，导数就代表一个变化率，从第一列变为第三列，灰度值相减，当然就是一个变化率了。这就是所谓的Prewitt算子。这样近似X方向导数就求出来了。Y方向导数与X方向导数求法相似，只不过是用第三行元素之和减去第一行元素之和。X方向和Y方向导数有了，那么梯度也就出来了。这样就可以找出一幅图中的边缘了。

还有一个问题，由于求的是3×3中心点的导数，所以给第二列加了一个权重，它的权重为2，第一列和第三列的权重为1，好了，这就是Sobel算子了。相比Prewitt算子，Sobel的抗噪能力更强。”

关键词：
一阶微分算子，Robert, Prewitt, Sobel, Krisch, Isotropic Sobel（事实的真相：理论中的一阶微分求导在应用中并不能实现，所以将其近似转化为一个小的像素模版对图像扫描，边扫描边刷新像素值），边缘检测，锐化，图像增强，梯度，卷积（MD，深入解释要走题，下一篇再解释这个词）

亲，看到这里，是不是感觉所谓itk工具包中以外国人命名的若干算法实现并没有那么难？逐渐抛离工具包中的方法，只用了其定义的基本结构。

静下心来，这些都属于玄门中的“术”，术业有专攻，术而已，唯手熟尔~

参考文献：
1.http://www.vtk.org/Wiki/ITK/Examples/EdgesAndGradients/SobelEdgeDetectionImageFilter
2.http://www.vtk.org/Wiki/ITK/Examples/Operators/SobelOperator
3.http://blog.csdn.NET/goodshot/article/details/10169619
4.http://blog.sciencenet.cn/blog-425437-776050.html
5.http://blog.csdn.net/zkp0601/article/details/42709669
6.http://blog.csdn.net/TonyShengTan/article/details/43371031
7.http://blog.sina.com.cn/s/blog_82a927880102vd9p.html