itk中的特征提取算法（一）

在计算机视觉中，特征提取是一个很重要的领域。而图像中被提取目标的边缘是它的几个重要特征之一，边缘标定了目标的位置范围，经常被用作图像匹配。
下面是几种通用的边缘检测算子（直接copy来的，若有差错，请指出）：
一阶：Roberts（Cross）算子，Prewitt算子，Sobel算子， Kirsch算子，罗盘算子；

二阶：Marr-Hildreth（LOG）算子，Canny算子，Laplacian算子。

好了，这里不讲大道理，理论知识网上都说的很清楚了，主要和itk官方大牛学算法代码。本文带大家了解下itk中的边缘检测算法。

1.Canny算子

//---------canny ，最重要的三个参数
float variance = 1.0;
float maximumError = 0.10;
float upperThreshold = 0.5;
typedef itk::MedianImageFilter<ImageType, ImageType > FilterType;
FilterType::Pointer medianFilter = FilterType::New();
FilterType::InputSizeType radius;
radius.Fill(3);
medianFilter->SetRadius(radius);
medianFilter->SetInput( input_data);
//------------------------------------------------------------
typedef double        RealPixelType;  //  Operations 注意这里需要转换成double
typedef itk::Image< RealPixelType, ImageDimension> RealImageType;
typedef itk::CastImageFilter< ImageType, RealImageType >    CastToRealFilterType;    
typedef itk::CannyEdgeDetectionImageFilter< RealImageType, RealImageType > CannyFilterType;
typedef itk::RescaleIntensityImageFilter< RealImageType, ImageType > RescaleFilterType;
CastToRealFilterType::Pointer toReal      = CastToRealFilterType::New();
CannyFilterType::Pointer      cannyFilter = CannyFilterType::New();
RescaleFilterType::Pointer    rescale     = RescaleFilterType::New();
toReal->SetInput( medianFilter->GetOutput() );
cannyFilter->SetInput( toReal->GetOutput() );
rescale->SetInput( cannyFilter->GetOutput() );
cannyFilter->SetVariance( variance );
cannyFilter->SetUpperThreshold( upperThreshold );
cannyFilter->SetMaximumError( maximumError );
//------------------------------------------------------------
output_data = rescale->GetOutput();

2.Sobel算子

typedef itk::SobelEdgeDetectionImageFilter <ImageType, FloatImageType> ImageFilterType; 
ImageFilterType::Pointer sobelFilter  = ImageFilterType::New(); 
sobelFilter->SetInput(input_data);
output_data = sobelFilter->GetOutput(); 

This filter uses the Sobel operator to calculate the image gradient and then finds the magnitude of this gradient vector.  The Sobel gradient magnitude (square-root sum of squares) is an indication of edge strength.
接下来看下实现方法：

typedef NeighborhoodOperatorImageFilter<InputImageType, OutputImageType>  OpFilter; 
typedef MultiplyImageFilter<OutputImageType,  OutputImageType, OutputImageType> MultFilter; 
typedef NaryAddImageFilter<OutputImageType, OutputImageType> AddFilter; //累计求和方案
typedef SqrtImageFilter<OutputImageType, OutputImageType>    SqrtFilter; //
 
unsigned int i;  
 
typename TOutputImage::Pointer output = this->GetOutput();  
output->SetBufferedRegion(output->GetRequestedRegion());
output->Allocate();
 
// Create the sobel operator
SobelOperator<OutputPixelType, ImageDimension> opers[ImageDimension];
ZeroFluxNeumannBoundaryCondition<TOutputImage> nbc;
 
// Setup mini-pipelines along each axis.
typename OpFilter::Pointer opFilter[ImageDimension]; //这个写法比较帅啊~
typename MultFilter::Pointer multFilter[ImageDimension];
typename AddFilter::Pointer addFilter = AddFilter::New();
typename SqrtFilter::Pointer sqrtFilter =  SqrtFilter::New();  
for(i=0; i < ImageDimension; ++i)
{
    // Create the filters for this axis.
    opFilter[i] = OpFilter::New();
    multFilter[i] = MultFilter::New();
   
    // Set boundary condition and operator for this axis.
    opers[i].SetDirection(i);
    opers[i].CreateDirectional();
    opFilter[i]->OverrideBoundaryCondition(&nbc);
    opFilter[i]->SetOperator(opers[i]);
   
    // Setup the mini-pipeline for this axis.
    opFilter[i]->SetInput(this->GetInput());
    multFilter[i]->SetInput1(opFilter[i]->GetOutput()); //这个平方。。。
    multFilter[i]->SetInput2(opFilter[i]->GetOutput());
   
    // All axes' mini-pipelines come together in addFilter.
    addFilter->SetInput(i, multFilter[i]->GetOutput());
} 
// calculate the gradient magnitude
sqrtFilter->SetInput(addFilter->GetOutput()); //开方
 
// setup the mini-pipeline to calculate the correct regions and
// write to the appropriate bulk data block
sqrtFilter->GraftOutput( this->GetOutput() );
 
// execute the mini-pipeline
sqrtFilter->Update();

3.Zero-crossing算子

double var = 5.0;
typedef itk::ZeroCrossingBasedEdgeDetectionImageFilter<ImageType, FloatImageType >FilterType; 
FilterType::Pointer edgeDetector = FilterType::New(); 
edgeDetector->SetInput( reader->GetOutput() ); 
  
FilterType::ArrayType variance;
variance.Fill(var); 
 
edgeDetector->SetVariance(variance);

ZeroCrossingBasedEdgeDetectionImageFilter方法中有个最重要的头文件：itkZeroCrossingImageFilter，提取出零交叉点。


“这个算子设计的初衷是用在那些通过计算二阶导数的边缘检测算子的后面，因为在计算出二阶导数之后，再用zero-crossing检测出来的zero-crossings就是二阶导数为0的点，也就是边缘。在数学意义上讲，这种做法就是认为二阶导数为0的点为边缘点。 ”

自我感觉，边缘检测在2D医学图像中有一定应用范围，而在3D医学图像中并不是很好用，边缘检测依赖于阈值的取舍，3D医学图像中病灶形状的不确定性增加了边缘检测的难度。边缘检测和区域分割是图像分割的两种不同思路，怎么选择，看临床应用。我用的版本是2.4，现成的边缘检测方法就上面三个，其实把算法原理弄清楚，自己可以随意写任何一种，“君子生非异也，善假于物也~”，工具包只是工具，不要拘泥在其中，不止要会用现成的类，还要在其基础结构上自由发挥。

参考文献：

1.http://blog.csdn.net/xiaowei_cqu/article/details/7829481

2.https://en.wikipedia.org/wiki/Canny_edge_detector

3.http://blog.csdn.net/poem_qianmo/article/details/25560901

4.http://blog.sina.com.cn/s/blog_710e52270102uzoy.html