ITK 实例4 OTSU算法对PNG图像进行多阈值二维分割

发布时间 2023-08-16 14:52:52作者: 一杯清酒邀明月
  1 #include "itkOtsuMultipleThresholdsCalculator.h"//包含头文件
  2  
  3 #include "itkImage.h"
  4 #include "itkImageFileReader.h"
  5 #include "itkImageFileWriter.h"
  6 #include "itkScalarImageToHistogramGenerator.h"
  7 #include "itkBinaryThresholdImageFilter.h"
  8 #include "itkNumericTraits.h"
  9  
 10 #include <iomanip>
 11 #include <stdio.h>
 12 //这个例子阐述了如何使用 itk::OtsuMultipleThresholdsCalculator
 13 int main( int argc, char * argv[] )
 14 {
 15   /*if( argc < 5 )
 16     {
 17     std::cerr << "Usage: " << argv[0];
 18     std::cerr << " inputImageFile outputImageFileBase ";
 19     std::cerr << "  outputImageFileExtension numberOfThresholdsToCalculate "  << std::endl;
 20     return EXIT_FAILURE;
 21     }*/
 22  
 23   //Convenience typedefs
 24   typedef  unsigned short  InputPixelType;
 25   typedef  unsigned char   OutputPixelType;
 26  
 27   typedef itk::Image< InputPixelType,  2 >   InputImageType;
 28   typedef itk::Image< OutputPixelType, 2 >   OutputImageType;
 29  
 30   /*OtsuMultipleThresholdsCalculator 以最大化类之间的变异的原则来计算一个给定的直方
 31   图的门限。我们使用 ScalarImageToHistogramGenerator 来得到直方图*/
 32   typedef itk::Statistics::ScalarImageToHistogramGenerator<
 33    InputImageType > ScalarImageToHistogramGeneratorType;
 34  
 35   typedef ScalarImageToHistogramGeneratorType::HistogramType HistogramType;
 36  
 37   typedef itk::OtsuMultipleThresholdsCalculator< HistogramType > CalculatorType;
 38   
 39   typedef itk::ImageFileReader< InputImageType >  ReaderType;
 40   typedef itk::ImageFileWriter< OutputImageType > WriterType;
 41   //一旦计算得到门限值,我们将使用 BinaryThresholdImageFilter 来对图像进行分割
 42    typedef itk::BinaryThresholdImageFilter<
 43   InputImageType, OutputImageType >  FilterType;
 44   
 45   ScalarImageToHistogramGeneratorType::Pointer scalarImageToHistogramGenerator
 46     = ScalarImageToHistogramGeneratorType::New();
 47  
 48   CalculatorType::Pointer calculator = CalculatorType::New();
 49   FilterType::Pointer filter = FilterType::New();
 50  
 51   ReaderType::Pointer reader = ReaderType::New();
 52   WriterType::Pointer writer = WriterType::New();
 53  
 54   scalarImageToHistogramGenerator->SetNumberOfBins( 128 );
 55   //设置Otsu分割阈值数目
 56   calculator->SetNumberOfThresholds( atoi("3") );
 57  
 58   const OutputPixelType outsideValue = 0;
 59   const OutputPixelType insideValue = 255;
 60  
 61   filter->SetOutsideValue( outsideValue );
 62   filter->SetInsideValue(  insideValue  );
 63  
 64   //Connect Pipeline
 65   reader->SetFileName( "BrainProtonDensitySlice.png" );
 66  
 67   //管道如下所示
 68   scalarImageToHistogramGenerator->SetInput( reader->GetOutput() );
 69   calculator->SetInputHistogram(scalarImageToHistogramGenerator->GetOutput() );
 70   filter->SetInput( reader->GetOutput() );
 71   writer->SetInput( filter->GetOutput() );
 72   //读取报错异常
 73   try
 74     {
 75     reader->Update();
 76     }
 77   catch( itk::ExceptionObject & excp )
 78     {
 79     std::cerr << "Exception thrown while reading image" << excp << std::endl;
 80     }
 81   scalarImageToHistogramGenerator->Compute();
 82   //calculator报错异常
 83   try
 84     {
 85     calculator->Compute();
 86     }
 87   catch( itk::ExceptionObject & excp )
 88     {
 89     std::cerr << "Exception thrown " << excp << std::endl;
 90     }
 91  
 92   //使用 GetOutput 方式可以得到门限值
 93   const CalculatorType::OutputType &thresholdVector = calculator->GetOutput();
 94  
 95   std::string outputFileBase = "BrainProtonDensitySlice_OtsuMu";
 96  
 97   InputPixelType lowerThreshold = itk::NumericTraits<InputPixelType>::min();
 98   InputPixelType upperThreshold;
 99  
100   typedef CalculatorType::OutputType::const_iterator ThresholdItType;
101  
102   for( ThresholdItType itNum = thresholdVector.begin();
103        itNum != thresholdVector.end();
104        ++itNum )
105     {
106     std::cout << "OtsuThreshold["
107               << (int)(itNum - thresholdVector.begin())
108               << "] = "
109               << static_cast<itk::NumericTraits<
110                           CalculatorType::MeasurementType>::PrintType>(*itNum)
111               << std::endl;
112  
113     upperThreshold = static_cast<InputPixelType>(*itNum);
114  
115     filter->SetLowerThreshold( lowerThreshold );
116     filter->SetUpperThreshold( upperThreshold );
117  
118     lowerThreshold = upperThreshold;
119  
120     std::ostringstream outputFilename;
121     outputFilename << outputFileBase
122                    << std::setfill('0') << std::setw(3) << (itNum - thresholdVector.begin())
123                    << "."
124                    << "png";
125     writer->SetFileName( outputFilename.str() );
126  
127     try
128       {
129       writer->Update();
130       }
131     catch( itk::ExceptionObject & excp )
132       {
133       std::cerr << "Exception thrown " << excp << std::endl;
134       }
135     }
136  
137   upperThreshold = itk::NumericTraits<InputPixelType>::max();
138   filter->SetLowerThreshold( lowerThreshold );
139   filter->SetUpperThreshold( upperThreshold );
140  
141   std::ostringstream outputFilename2;
142   outputFilename2 << outputFileBase
143                   << std::setfill('0') << std::setw(3) << thresholdVector.size()
144                   << "."
145                   << "jpg";
146   writer->SetFileName( outputFilename2.str() );
147  
148   try
149     {
150     writer->Update();
151     }
152   catch( itk::ExceptionObject & excp )
153     {
154     std::cerr << "Exception thrown " << excp << std::endl;
155     }
156  
157   return EXIT_SUCCESS;
158 }