这篇算是《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,侵删):
关键词:
一阶微分算子,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