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

本文的重点是图像处理中经典hough变换之直线检测，类名叫做itkHoughTransform2DLinesImageFilter（它还有个兄弟叫做itkHoughTransform2DCirclesImageFilter，下次再讲）。

描述：执行Hough变换在二维图像中找到二维直线。
该过滤器源自ImageToImageFilter；
输入：一个图像，通过阈值提取需要的像素点；
输出1：累加器构成的图像；（注意这里，不是输出一幅画好直线的图像，并不是直接输出结果）
输出2：GetLines()：返回了一个LinesSpatialObjects的列表；（这个输出也不是图像，而是直线的各种相关信息）
形状参数化后的直线方程：R = x*vcl_cos(Teta)+y*vcl_sin(Teta)；（R是到原点的垂直距离，TETA是法线的夹角）

今天的使用方法有点复杂~：

	typedef   double                                   HoughSpacePixelType;
	typedef   itk::Image< HoughSpacePixelType, 2>      HoughImageType;

	typedef itk::CastImageFilter<ImageType, HoughImageType >    CastingFilterType;
	CastingFilterType::Pointer caster = CastingFilterType::New();
	caster->SetInput(reader->GetOutput());

	typedef itk::GradientMagnitudeImageFilter<HoughImageType,HoughImageType> GradientFilterType;
	GradientFilterType::Pointer gradFilter =  GradientFilterType::New();
	gradFilter->SetInput(caster->GetOutput());
	gradFilter->Update();

	typedef itk::ThresholdImageFilter<HoughImageType> ThresholdFilterType;
	ThresholdFilterType::Pointer threshFilter = ThresholdFilterType::New();
	threshFilter->SetInput(gradFilter->GetOutput());
	threshFilter->SetOutsideValue(0);
	unsigned char thresh_below = 10;
	unsigned char thresh_above = 200;
	threshFilter->ThresholdOutside(thresh_below,thresh_above);
	threshFilter->Update();

	typedef itk::HoughTransform2DLinesImageFilter<HoughSpacePixelType,HoughSpacePixelType> HoughTransformFilterType;  
	HoughTransformFilterType::Pointer houghFilter = HoughTransformFilterType::New();
	houghFilter->SetInput(threshFilter->GetOutput());  
	houghFilter->SetThreshold(0.0f);
	houghFilter->SetAngleResolution(500.0f);
	houghFilter->SetDiscRadius(10.0f);
	houghFilter->SetVariance(10.0f);
	houghFilter->Update();
	houghFilter->Simplify();

接下来源码分析：
1..h中的私有变量：

  float              m_AngleResolution;//角分辨：就是这里 500
  float              m_Threshold;//阈值？ 0
  OutputImagePointer m_SimplifyAccumulator;//简单的迭代器 NULL
  LinesListType      m_LinesList;//线的空间对象
  unsigned int       m_NumberOfLines;//线的条数 1
  float              m_DiscRadius;//半径 10
  float              m_Variance;//方差 5
  unsigned long      m_OldModifiedTime;//时间 0
  unsigned long      m_OldNumberOfLines;//线的条数？ 0

2.生成输出图像详细信息的函数：

template<typename TInputPixelType, typename TOutputPixelType>
void 
HoughTransform2DLinesImageFilter<TInputPixelType,TOutputPixelType>
::GenerateOutputInformation()
{
  // call the superclass' implementation of this method
  Superclass::GenerateOutputInformation();

  // get pointers to the input and output
  InputImageConstPointer  input  = this->GetInput();
  OutputImagePointer      output = this->GetOutput();

  if ( !input || !output )
    {
    return;
    }

  // 注意这里Compute the size of the output image
  typename InputImageType::RegionType region;
  Size<2> size;
  size[0]= (long unsigned int)(vcl_sqrt(m_AngleResolution*m_AngleResolution+input->GetLargestPossibleRegion().GetSize()[0]*input->GetLargestPossibleRegion().GetSize()[0]));
  size[1]= (long unsigned int)m_AngleResolution;
  region.SetSize(size);
  region.SetIndex(input->GetLargestPossibleRegion().GetIndex());

  output->SetLargestPossibleRegion( region );
}

3.数据生成

template<typename TInputPixelType, typename TOutputPixelType>
void
HoughTransform2DLinesImageFilter< TInputPixelType, TOutputPixelType>
::GenerateData()
{
  itkDebugMacro(<<"HoughTransform2DLinesImageFilter called");

  // Get the input and output pointers
  InputImageConstPointer  inputImage  = this->GetInput(0);
  OutputImagePointer outputImage = this->GetOutput(0);

  // Allocate the output
  this->AllocateOutputs();
  outputImage->FillBuffer(0);

  const double nPI = 4.0 * vcl_atan( 1.0 );
//又见像素迭代器，几乎每个算法类都能见到
  ImageRegionConstIteratorWithIndex< InputImageType >  image_it( inputImage,  inputImage->GetRequestedRegion() );
  image_it.Begin();

  Index<2> index;

  while( !image_it.IsAtEnd() )
    {
    if(image_it.Get()>m_Threshold)
      { 
      for(double angle = -nPI; angle < nPI; angle += nPI / m_AngleResolution )
        {  
//公式在这里
        index[0]=(long unsigned int)(image_it.GetIndex()[0]*vcl_cos(angle)+image_it.GetIndex()[1]*vcl_sin(angle)); // m_R
        index[1]= (long unsigned int)((m_AngleResolution/2)+m_AngleResolution*angle/(2*nPI)); // m_Theta
//在符合条件的点上做标记  
        if ( (index[0] > 0) &&
             (index[0] <= (long)outputImage->GetBufferedRegion().GetSize()[0]))
          // the preceeding "if" should be replacable with "if (
          // outputImage->GetBufferedRegion().IsInside(index) )" but
          // the algorithm fails if it is 
          {
          outputImage->SetPixel(index, outputImage->GetPixel(index)+1);  
          }
        } 
      }
    ++image_it;
    }
}

4. 头文件中include 的itkImageToImageFilter
一直好奇这个类有什么用，现在看起来简直堪称神作啊~
算法类都是这个类的子类。子类中只重写GenerateData就差不多了。

最近心太浮躁了，去看看C++面向对象。

参考文献：
http://blog.csdn.net/abcjennifer/article/details/7448513
http://blog.csdn.net/carson2005/article/details/6568414
http://blog.163.com/yuyang_tech/blog/static/21605008320130233343990/
http://blog.sina.com.cn/s/blog_500bd63c0101fdmx.html
http://www.docin.com/p-67248749.html