int
main(int argc, char * argv[])
{
if (argc < 3)
{
std::cerr << "Usage: " << std::endl;
std::cerr << argv[0];
std::cerr << " <InputFileName> <OutputFileName>";
std::cerr << " [Sigma]";
std::cerr << std::endl;
return EXIT_FAILURE;
}
const char * inputFileName = argv[1];
const char * outputFileName = argv[2];
double sigmaValue = 0.5;
if (argc > 3)
{
sigmaValue = std::stod(argv[3]);
}
using PixelType = float;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(inputFileName);
PadFilterType::Pointer padFilter = PadFilterType::New();
padFilter->SetInput(reader->GetOutput());
padding[0] = 0;
padding[1] = 2;
padding[2] = 6;
padFilter->SetPadUpperBound(padding);
using ComplexImageType = ForwardFFTFilterType::OutputImageType;
ForwardFFTFilterType::Pointer forwardFFTFilter = ForwardFFTFilterType::New();
forwardFFTFilter->SetInput(padFilter->GetOutput());
try
{
forwardFFTFilter->UpdateOutputInformation();
}
catch (itk::ExceptionObject & error)
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
GaussianSourceType::Pointer gaussianSource = GaussianSourceType::New();
gaussianSource->SetNormalized(true);
ComplexImageType::ConstPointer transformedInput = forwardFFTFilter->GetOutput();
const ComplexImageType::SpacingType inputSpacing = transformedInput->GetSpacing();
gaussianSource->SetSize(inputSize);
gaussianSource->SetSpacing(inputSpacing);
gaussianSource->SetOrigin(inputOrigin);
gaussianSource->SetDirection(inputDirection);
GaussianSourceType::ArrayType sigma;
for (
unsigned int ii = 0; ii <
Dimension; ++ii)
{
const double halfLength = inputSize[ii] * inputSpacing[ii] / 2.0;
sigma[ii] *= halfLength;
mean[ii] = inputOrigin[ii] + halfLength;
}
mean = inputDirection * mean;
gaussianSource->SetSigma(sigma);
gaussianSource->SetMean(mean);
FFTShiftFilterType::Pointer fftShiftFilter = FFTShiftFilterType::New();
fftShiftFilter->SetInput(gaussianSource->GetOutput());
MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New();
multiplyFilter->SetInput1(forwardFFTFilter->GetOutput());
multiplyFilter->SetInput2(fftShiftFilter->GetOutput());
InverseFilterType::Pointer inverseFFTFilter = InverseFilterType::New();
inverseFFTFilter->SetInput(multiplyFilter->GetOutput());
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(outputFileName);
writer->SetInput(inverseFFTFilter->GetOutput());
try
{
writer->Update();
}
catch (itk::ExceptionObject & error)
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}