[Insight-users] Adaptive thresholding by region

[Loïc Giraldi]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100510/b6757726/attachment-0001.htm>