itkRichardsonLucyDeconv.cxx

#include "itkFFTConvolutionImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRichardsonLucyDeconvolutionImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkMacro.h"
#include "itkMath.h"

int main(int argc, char* argv[])
{
    if (argc < 6) {
        std::cerr << "Usage: " << argv[0]
            << " <input image> <kernel image> <blur image> <output image> <iterations> <relaxation>"
            << std::endl;
        return EXIT_FAILURE;
    }

    typedef float                                               RealPixelType;
    typedef unsigned char                                       UCharPixelType;
    const unsigned int                                          Dimension = 2;
    typedef itk::Image< RealPixelType, Dimension >              ImageType;
    typedef itk::Image<UCharPixelType, Dimension >              UCharImageType;
    typedef itk::ImageFileReader< ImageType >                   ReaderType;
    typedef itk::ImageFileWriter<UCharImageType >               WriterType;
    typedef itk::FFTConvolutionImageFilter< ImageType >         ConvolutionFilterType;
    typedef itk::RichardsonLucyDeconvolutionImageFilter< ImageType > DeconvolutionFilterType;
    typedef itk::RescaleIntensityImageFilter< ImageType, UCharImageType > RescaleFilterType;

    ReaderType::Pointer inputReader = ReaderType::New();
    inputReader->SetFileName(argv[1]);
    inputReader->Update();

    ReaderType::Pointer kernelReader = ReaderType::New();
    kernelReader->SetFileName(argv[2]);
    kernelReader->Update();

    // Generate a convolution of the input image with the kernel image
    ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New();
    convolutionFilter->SetInput(inputReader->GetOutput());
    convolutionFilter->NormalizeOn();
    convolutionFilter->SetKernelImage(kernelReader->GetOutput());

    // Test the deconvolution algorithm
    DeconvolutionFilterType::Pointer deconvolutionFilter = DeconvolutionFilterType::New();
    deconvolutionFilter->SetInput(convolutionFilter->GetOutput());
    deconvolutionFilter->SetKernelImage(kernelReader->GetOutput());
    deconvolutionFilter->NormalizeOn();

    unsigned int iterations = static_cast< unsigned int >(atoi(argv[5]));
    deconvolutionFilter->SetNumberOfIterations(iterations);

    RescaleFilterType::Pointer realRescaleInputFilter = RescaleFilterType::New();
    realRescaleInputFilter->SetInput(convolutionFilter->GetOutput());
    realRescaleInputFilter->SetOutputMinimum( itk::NumericTraits< UCharPixelType >::min() );
    realRescaleInputFilter->SetOutputMaximum( itk::NumericTraits< UCharPixelType >::max() );

    RescaleFilterType::Pointer realRescaleOutputFilter = RescaleFilterType::New();
    realRescaleOutputFilter->SetInput(deconvolutionFilter->GetOutput());
    realRescaleOutputFilter->SetOutputMinimum( itk::NumericTraits< UCharPixelType >::min() );
    realRescaleOutputFilter->SetOutputMaximum( itk::NumericTraits< UCharPixelType >::max() );

    // Write the deconvolution result
    try {
        WriterType::Pointer writer = WriterType::New();

        writer->SetFileName(argv[3]);
        writer->SetInput(realRescaleInputFilter->GetOutput());
        writer->Update();

        writer->SetFileName(argv[4] );
        writer->SetInput(realRescaleOutputFilter->GetOutput());
        writer->Update();
    } catch ( itk::ExceptionObject & e ) {
        std::cerr << "Unexpected exception caught when writing deconvolution image: "
            << e << std::endl;
        return EXIT_FAILURE;
    }

    return EXIT_SUCCESS;
}