itkTikhonov.cxx

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

int main(int argc, char* argv[])
{
    if (argc < 5) {
        std::cerr << "Usage: " << argv[0]
            << " <input image> <kernel image> <blur image> <output image> <noise_variance>"
            << 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::TikhonovDeconvolutionImageFilter< 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->SetNormalize(true);

    double noise_variance = static_cast<double>(atof(argv[5]));
    deconvolutionFilter->SetRegularizationConstant(noise_variance);

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

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

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

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

        writer->SetFileName(argv[4] );
        writer->SetInput(rescaleOutputFilter->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;
}