[Insight-users] Adaptive thresholding by region

[Lassi Paavolainen]

Hi Loic,

You might want to have a look at following Insight Journal paper that 
should do exactly what you plan to do:
http://insight-journal.org/browse/publication/702

Best regards,
Lassi

On Mon, 10 May 2010, Loïc Giraldi wrote:

> Hi all!
> 
> I'm trying to make an adaptive thresholding using the Otsu method. To do so,
> I want to process my image by region. The idea is to separate the original
> picture in several blocks, process each block and reconstruct the final
> output by gathering the blocks.
> I am trying to do it using the RegionOfInterestImageFilter, the
> OtsuThresholdImageFilter and the PasteImageFilter in a loop over regions, I
> have grabbed some informations in this mailing list but I can't manage to
> make it working. The process works if I do not loop over regions.
> 
> Here is my code:
> 
> #include "itkImage.h"
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
> 
> #include "itkRegionOfInterestImageFilter.h"
> #include "itkOtsuThresholdImageFilter.h"
> #include "itkPasteImageFilter.h"
> 
> #include <iostream>
> #include <string>
> #include <sstream>
> #include <cstdlib>
> using namespace std;
> 
> 
> template <typename T> string toString(T val){
>     ostringstream oss;
>     oss << val;
>     return oss.str();
> }
> 
> int main(int argc, char *argv[]){
>     ////////////////////////////  Usage //////////////////////////////
>     if (argc<3){
>         std::cerr << "Missing arguments!" << std::endl;
>         std::cerr << "Usage : " << argv[0]
>                   << " InputFile"
>                   << " Size" << std::endl;
>         return 1;
>     }
>        
>     const char *inputFileName = argv[1];
>     const unsigned int size = atoi(argv[2]);
>    
>     //////////////////////  Image definitions ////////////////////////
>     const int Dimension = 3;
> 
>     typedef unsigned short InputPixelType;
>     typedef unsigned char OutputPixelType;
> 
>     typedef itk::Image< InputPixelType, Dimension > InputImageType;
>     typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
> 
>     ////////////////////////  Image reading //////////////////////////
>     typedef itk::ImageFileReader< InputImageType > ReaderType;
>    
>     ReaderType::Pointer reader = ReaderType::New();
>     reader->SetFileName(inputFileName);
>     try{
>         reader->Update();
>     }
>     catch( itk::ExceptionObject & excep ){
>         std::cerr << "Exception caught !" << std::endl;
>         std::cerr << excep << std::endl;
>     }
>    
>     //////////////////////////  Filtering ////////////////////////////
>     //Output Allocation
>     OutputImageType::Pointer output = OutputImageType::New();
>     OutputImageType::IndexType outIndex;
>     OutputImageType::SizeType outSize;
>     OutputImageType::RegionType outRegion;
> 
>     for (int i = 0; i < Dimension; ++i){
>         outIndex[i]=reader->GetOutput()->GetBufferedRegion().GetIndex()[i];
>         outSize[i]=reader->GetOutput()->GetBufferedRegion().GetSize()[i];
>     }
> 
>     outRegion.SetIndex(outIndex);
>     outRegion.SetSize(outSize);
>    
>     output->SetRegions(outRegion);
>     output->Allocate();
>     output->FillBuffer(0);
>    
>     //Moving box definition
>     InputImageType::IndexType boxIndex;
>     InputImageType::SizeType boxSize;
>     InputImageType::RegionType boxRegion;
> 
>     for (int i = 0; i < Dimension; ++i){
>         boxIndex[i] = outIndex[i];
>     }
>     for (int i = 0; i < Dimension; ++i){
>         boxSize[i] = size;
>     }
>    
>     boxRegion.SetIndex(boxIndex);
>     boxRegion.SetSize(boxSize);
> 
>     unsigned numberOfBoxes[3];
>     for (int i = 0; i < Dimension; ++i){
>         numberOfBoxes[i] = static_cast<unsigned>((outSize[i]/size));
>     }
> 
>     //Processing
>     typedef itk::RegionOfInterestImageFilter< InputImageType,
> InputImageType> ROIFilterType;
>     typedef itk::OtsuThresholdImageFilter< InputImageType, OutputImageType >
> OtsuThresholdFilterType;
>     typedef itk::PasteImageFilter<OutputImageType> PasteImageFilterType;
>    
>     ROIFilterType::Pointer roi = ROIFilterType::New();
>     OtsuThresholdFilterType::Pointer otsu = OtsuThresholdFilterType::New();
>     PasteImageFilterType::Pointer paste = PasteImageFilterType::New();
> 
>     for(unsigned i = 0; i < numberOfBoxes[0]; ++i){
>         for(unsigned j = 0; j < numberOfBoxes[1]; ++j){
>             for(unsigned k = 0; k < numberOfBoxes[2]; ++k){
>                
>                 boxIndex[0] += i * size;
>                 boxIndex[1] += j * size;
>                 boxIndex[2] += k * size;
>                 boxRegion.SetIndex(boxIndex);
>                
>                 roi->SetInput(reader->GetOutput());
>                 roi->SetRegionOfInterest(boxRegion);
> 
>                 try{
>                     roi->UpdateLargestPossibleRegion();
>                     roi->Update();
>                 }
>                 catch( itk::ExceptionObject & excep ){
>                     std::cerr << "ROI!" << std::endl;
>                     std::cerr << "Exception caught !" << std::endl;
>                     std::cerr << excep << std::endl;
>                 }
>                
>                 otsu->SetInput(roi->GetOutput());
>                 otsu->SetNumberOfHistogramBins(128);
>                 otsu->SetInsideValue(255);
>                 otsu->SetOutsideValue(0);
>                
>                 try{
>                     otsu->UpdateLargestPossibleRegion();
>                     otsu->Update();
>                 }
>                 catch( itk::ExceptionObject & excep ){
>                     std::cerr << "Otsu!" << std::endl;
>                     std::cerr << "Exception caught !" << std::endl;
>                     std::cerr << excep << std::endl;
>                 }
>                 otsu->Print(std::cout);
>                
>                 paste->SetDestinationImage(output);
>                 paste->SetSourceImage(otsu->GetOutput());
>                 paste->SetDestinationIndex( boxIndex );
>                 paste->SetSourceRegion( boxRegion );
> 
>                 try{
>                     paste->UpdateLargestPossibleRegion();
>                     paste->Update();
>                 }
>                 catch( itk::ExceptionObject & excep ){
>                     std::cerr << "Paste!" << std::endl;
>                     std::cerr << "Exception caught !" << std::endl;
>                     std::cerr << excep << std::endl;
>                     std::cerr << "i,j,k = " << i << ',' << j << ',' << k <<
> std::endl;
>                     boxRegion.Print(std::cerr);
>                     paste.Print(std::cerr);
>                     paste->GetOutput()->Print(std::cerr);
>                 }
>                 output = paste->GetOutput();
>                 output->DisconnectPipeline();
>             }
>         }
>     }
>    
>     ////////////////////////  Image writing //////////////////////////
>     typedef itk::ImageFileWriter< OutputImageType > WriterType;
> 
>     WriterType::Pointer writer = WriterType::New();
> 
>     string outputFileName(inputFileName);
>     size_t pos = outputFileName.find_last_of('.');
>     outputFileName.insert(pos, "_result_"+toString(size));
> 
>     writer->SetInput(output);
>     writer->SetFileName(outputFileName.c_str());
>     writer->Update();
>     try{
>         writer->Update();
>     }
>     catch( itk::ExceptionObject & excep ){
>         std::cerr << "Write!" << std::endl;
>         std::cerr << "Exception caught !" << std::endl;
>         std::cerr << excep << std::endl;
>     }
>     return 0;
> }
> 
> Thanks for your help,
> Regards,
> Loïc
> 
> 
>