ITK/Examples/WishList/Segmentation/OtsuMultipleThresholdsCalculator

/* Just a script to make one image out of otsu's multiple threshold scheme. Say, if one defines two threshold values, then there will be an image with three values: value 1 for (min value of the image < first threshold),( first threshold < value 2. < second threshold) and (second threshold < value 3 < max value of the image). Just change the imagetype for different images.

example:

OtsuMultipleThresholdCalculator.exe inputimagename.png outputimagename .png 128 2

• /

1. include "itkOtsuMultipleThresholdsCalculator.h"
2. include "itkImage.h"
3. include "itkImageFileReader.h"
4. include "itkImageFileWriter.h"
5. include "itkScalarImageToHistogramGenerator.h"
6. include "itkBinaryThresholdImageFilter.h"
7. include "itkNumericTraits.h"
8. include <itkAddImageFilter.h>
9. include <stdio.h>
10. include <math.h>

typedef unsigned char PixelType; typedef itk::Image<PixelType, 2> ImageType;

void CreateImage(ImageType::Pointer image);

int main( int argc, char * argv[]) {

  if( argc < 6)
   {
   std::cerr << "Usage: " << argv[0];
   std::cerr << " inputImageFile outputImageFileBase ";  
   std::cerr << "  outputImageFileExtension numberOfhistogrambins " ;

std::cerr << " numberOfThresholdsToCalculate "<< std::endl;

   return EXIT_FAILURE;
   }

 typedef itk::ImageFileReader< ImageType >  ReaderType;
 ReaderType::Pointer reader = ReaderType::New();
 reader->SetFileName( argv[1] );
 reader->Update();

 ImageType::Pointer image = reader -> GetOutput();

 ImageType::Pointer sumimage = ImageType::New();

 ImageType::RegionType region = image->GetLargestPossibleRegion();

 ImageType::SizeType size = region.GetSize();

 ImageType::IndexType start;
 start.Fill(0);

 region.SetIndex( start );
 region.SetSize( size );

 sumimage->SetRegions(region);
 sumimage->Allocate();
 sumimage->FillBuffer(0) ;

 typedef itk::ImageFileWriter< ImageType >  WriterType; 
 typedef itk::Statistics::ScalarImageToHistogramGenerator<ImageType >
 ScalarImageToHistogramGeneratorType;
 typedef
 itk::OtsuMultipleThresholdsCalculator<ScalarImageToHistogramGeneratorType::HistogramType
 > CalculatorType;
 typedef itk::BinaryThresholdImageFilter< ImageType, ImageType >
 FilterType;
 typedef itk::AddImageFilter< ImageType, ImageType ,ImageType>
 AddFilterType;

 ScalarImageToHistogramGeneratorType::Pointer
 scalarImageToHistogramGenerator =
 ScalarImageToHistogramGeneratorType::New();
 CalculatorType::Pointer calculator = CalculatorType::New();
 FilterType::Pointer filter = FilterType::New();
 AddFilterType::Pointer addFilter = AddFilterType::New();
 
 

 WriterType::Pointer writer = WriterType::New();

 scalarImageToHistogramGenerator->SetNumberOfBins( atoi( argv[4] ) );
 calculator->SetNumberOfThresholds(atoi( argv[5] ) );

 scalarImageToHistogramGenerator->SetInput(image);
 scalarImageToHistogramGenerator->Compute();
 calculator->SetInputHistogram(scalarImageToHistogramGenerator->GetOutput());
 filter->SetInput(image );

 calculator->Update();

 const CalculatorType::OutputType &thresholdVector =
 calculator->GetOutput();
 CalculatorType::OutputType::const_iterator itNum =
 thresholdVector.begin();

 PixelType min = 0;
 PixelType max = 255;

 const PixelType outsideValue = 0;

 std::string outputFileBase = argv[2];
 std::string outputFile;
 std::string format = argv[2];

 char outputFilename[1000];
 outputFile = outputFileBase;
 outputFile += argv[3];   // filename extension

 PixelType lowerThreshold = min;
 PixelType upperThreshold;

 filter->SetOutsideValue(outsideValue );
 PixelType step =floor(255/(atoi( argv[5] )+1) );
 PixelType count=step;

 for(; itNum < thresholdVector.end(); itNum++)
 {

 const PixelType insideValue = count;
 filter->SetInsideValue( insideValue );

 upperThreshold =static_cast<PixelType>(*itNum);
 
 filter->SetLowerThreshold( lowerThreshold );
 filter->SetUpperThreshold(upperThreshold );
 filter->Update();

 lowerThreshold = upperThreshold;

 addFilter->SetInput1( filter->GetOutput() );
 addFilter->SetInput2(  sumimage );

 addFilter->Update();
 sumimage = addFilter->GetOutput() ;


 count = count + step;
 }

 PixelType insideValue = 255;


 filter->SetInsideValue( insideValue   );
 filter->SetLowerThreshold( lowerThreshold );
 filter->SetUpperThreshold(max );
 filter->Update();


 addFilter->SetInput1( filter->GetOutput() );
 addFilter->SetInput2(  sumimage );

 addFilter->Update();
 
 sumimage = addFilter->GetOutput() ;
 writer->SetInput(sumimage );

 writer->SetFileName(  outputFile.c_str());
   
   try
     { 
     writer->Update(); 
     }
   catch( itk::ExceptionObject & excp )
     {
     std::cerr << "Exception thrown " << excp << std::endl;
     }
   // Software Guide : BeginCodeSnippet


 return EXIT_SUCCESS;

}

CMakeLists.txt

cmake_minimum_required(VERSION 2.6)

project(OtsuMultipleThresholdsCalculator)

find_package(ITK REQUIRED) include(${ITK_USE_FILE})

add_executable(OtsuMultipleThresholdsCalculator OtsuMultipleThresholdsCalculator.cxx) TARGET_LINK_LIBRARIES(OtsuMultipleThresholdsCalculator ITKIO ITKStatistics)