itk中的图像分割算法（二）

今天要介绍的是经典的自动化阈值分割算法--Otsu，中文叫做大津算法（哎，就是日本的大津发明的这个算法，最烦这些翻译过来的鬼名字，云里雾里，直接叫最大类间方差法多好~）。

itk中提供的类：itkOtsuThresholdImageFilter

用法：

typedef itk::OtsuThresholdImageFilter<ImageType, ImageType >  FilterType;filter->SetInput( input_data );filter->SetNumberOfHistogramBins (atoi(argv[3]));//这个参数很重要filter->Update();std::cout << "Computed Threshold is: " << filter->GetThreshold() << std::endl;output_data = filter->GetOutput();

接下来带大家看源码：

0.include头文件

#include "itkBinaryThresholdImageFilter.h"//otsu的本质就是二值化图像#include "itkOtsuThresholdImageCalculator.h"//这里有个鬼？？？

1.私有变量

  InputPixelType      m_Threshold;  OutputPixelType     m_InsideValue;  OutputPixelType     m_OutsideValue;  unsigned long       m_NumberOfHistogramBins;

2.数据生成方法

template<class TInputImage, class TOutputImage>voidOtsuThresholdImageFilter<TInputImage, TOutputImage>::GenerateData(){  typename ProgressAccumulator::Pointer progress = ProgressAccumulator::New();//这部分照旧不看  progress->SetMiniPipelineFilter(this);  // Compute the Otsu Threshold for the input image 往这儿看@_@!!关键部分在这儿//这个类的重点就是计算出一个阈值：otsu->GetThreshold();  typename OtsuThresholdImageCalculator<TInputImage>::Pointer otsu =    OtsuThresholdImageCalculator<TInputImage>::New();  otsu->SetImage (this->GetInput());  otsu->SetNumberOfHistogramBins (m_NumberOfHistogramBins);  otsu->Compute();  m_Threshold = otsu->GetThreshold();  typename BinaryThresholdImageFilter<TInputImage,TOutputImage>::Pointer threshold =    BinaryThresholdImageFilter<TInputImage,TOutputImage>::New();  progress->RegisterInternalFilter(threshold,.5f);  threshold->GraftOutput (this->GetOutput());  threshold->SetInput (this->GetInput());  threshold->SetLowerThreshold(NumericTraits<InputPixelType>::NonpositiveMin());  threshold->SetUpperThreshold(otsu->GetThreshold());  threshold->SetInsideValue (m_InsideValue);  threshold->SetOutsideValue (m_OutsideValue);  threshold->Update();  this->GraftOutput(threshold->GetOutput());}

3.真相，，原来itkOtsuThresholdImageFilter只是个画皮呀。。。

用法：

typedef itk::OtsuThresholdImageCalculator<ImageType>  CalculatorType;    CalculatorType::Pointer calculator = CalculatorType::New();    calculator->SetImage(image);    calculator->SetNumberOfHistogramBins( 64);    calculator->Compute();    std::cout << "calculator: " << calculator;    std::cout << "NumberOfHistogramBins: " << calculator->GetNumberOfHistogramBins();    short thresholdResult = calculator->GetThreshold();    std::cout << "The threshold intensity value is : " << thresholdResult << std::endl;

4.深入解析OtsuThresholdImageCalculator计算方法，，，有点长

/* * Compute the Otsu's threshold */template<class TInputImage>voidOtsuThresholdImageCalculator<TInputImage>::Compute(void){  unsigned int j;  if ( !m_Image ) { return; } //防爆代码  if( !m_RegionSetByUser )//是否人工输入计算区域    {    m_Region = m_Image->GetRequestedRegion();//false的话，就获取整个图像区域    }  double totalPixels = (double) m_Region.GetNumberOfPixels();//获取图像的像素个数  if ( totalPixels == 0 ) { return; }//防爆代码  // compute image max and min //获取图像像素灰度的最大值和最小值  typedef MinimumMaximumImageCalculator<TInputImage> RangeCalculator;//最大最小值计算器。。。  typename RangeCalculator::Pointer rangeCalculator = RangeCalculator::New();  rangeCalculator->SetImage( m_Image );  rangeCalculator->Compute();  PixelType imageMin = rangeCalculator->GetMinimum();//最小值  PixelType imageMax = rangeCalculator->GetMaximum();//最大值  if ( imageMin >= imageMax )//两值相等时，那个大于没什么用，不过老鸟都会这么写，安全性    {    m_Threshold = imageMin;//直接将最小值当作结果返回    return;    }  // create a histogram  创建灰度直方图  std::vector<double> relativeFrequency;  relativeFrequency.resize( m_NumberOfHistogramBins );  for ( j = 0; j < m_NumberOfHistogramBins; j++ )    {    relativeFrequency[j] = 0.0;    }  double binMultiplier = (double) m_NumberOfHistogramBins /    (double) ( imageMax - imageMin );  typedef ImageRegionConstIteratorWithIndex<TInputImage> Iterator;  Iterator iter( m_Image, m_Region );  while ( !iter.IsAtEnd() )    {    unsigned int binNumber;    PixelType value = iter.Get();    if ( value == imageMin )       {      binNumber = 0;      }    else      {      binNumber = (unsigned int) vcl_ceil((value - imageMin) * binMultiplier ) - 1;      if ( binNumber == m_NumberOfHistogramBins ) // in case of rounding errors        {        binNumber -= 1;        }      }    relativeFrequency[binNumber] += 1.0;    ++iter;    }   // normalize the frequencies 将频率归一化，或者标准化  double totalMean = 0.0;  for ( j = 0; j < m_NumberOfHistogramBins; j++ )    {    relativeFrequency[j] /= totalPixels;    totalMean += (j+1) * relativeFrequency[j];    }  // compute Otsu's threshold by maximizing the between-class  // variance  double freqLeft = relativeFrequency[0];  double meanLeft = 1.0;  double meanRight = ( totalMean - freqLeft ) / ( 1.0 - freqLeft );  double maxVarBetween = freqLeft * ( 1.0 - freqLeft ) *    vnl_math_sqr( meanLeft - meanRight );  int maxBinNumber = 0;  double freqLeftOld = freqLeft;  double meanLeftOld = meanLeft;  for ( j = 1; j < m_NumberOfHistogramBins; j++ )    {    freqLeft += relativeFrequency[j];    meanLeft = ( meanLeftOld * freqLeftOld +                  (j+1) * relativeFrequency[j] ) / freqLeft;    if (freqLeft == 1.0)      {        meanRight = 0.0;      }    else      {      meanRight = ( totalMean - meanLeft * freqLeft ) /         ( 1.0 - freqLeft );      }    double varBetween = freqLeft * ( 1.0 - freqLeft ) *      vnl_math_sqr( meanLeft - meanRight );       if ( varBetween > maxVarBetween )      {      maxVarBetween = varBetween;      maxBinNumber = j;      }    // cache old values    freqLeftOld = freqLeft;    meanLeftOld = meanLeft;     }   m_Threshold = static_cast<PixelType>( imageMin +                                         ( maxBinNumber + 1 ) / binMultiplier );}

5.自己实现一个otsu，其实很简单么，没多少内容。

int Otsu(ImageType2D *src)  {  ImageType2D::RegionType inputRegion = src->GetLargestPossibleRegion();ImageType2D::SizeType size = inputRegion.GetSize();    int height= size[1];      int width= size[0];      //histogram float histogram[256] = {0}; typedef itk::ImageRegionIteratorWithIndex<ImageType2D>   FieldIterator;FieldIterator fieldIter( src, src->GetLargestPossibleRegion());fieldIter.GoToBegin();ImageType2D::IndexType index;while( !fieldIter.IsAtEnd()  ){index=fieldIter.GetIndex();unsigned char p = src->GetPixel(index);histogram[p]++;++fieldIter;}    //normalize histogram & average pixel value     int size_ = height * width;      float u =0;    for(int i = 0; i < 256; i++)    {          histogram[i] = histogram[i] / size_;          u += i * histogram[i];  //整幅图像的平均灰度    }      int threshold;        float maxVariance=0;      float w0 = 0, avgValue  = 0;    for(int i = 0; i < 256; i++)     {          w0 += histogram[i];  //假设当前灰度i为阈值, 0~i 灰度像素所占整幅图像的比例即前景比例        avgValue  += i * histogram[i]; //avgValue/w0 = u0        float t = avgValue/w0 - u;  //t=u0-u        float variance = t * t * w0 /(1 - w0);          if(variance > maxVariance)         {              maxVariance = variance;              threshold = i;          }      }     return threshold;  } 

最近一直在忙，这是这个月的第一篇，惭愧。

如何才能战无不胜？？？

自己需要熟悉还很多，合理安排时间和精力，继续吧~！

参考文献：

1.http://blog.csdn.net/u011285477/article/details/52004513