今天要介绍的是经典的自动化阈值分割算法--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>
void
OtsuThresholdImageFilter<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>
void
OtsuThresholdImageCalculator<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