/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

#ifndef itkFFTConvolutionImageFilter_hxx
#define itkFFTConvolutionImageFilter_hxx

#include "itkCastImageFilter.h"
#include "itkChangeInformationImageFilter.h"
#include "itkConstantPadImageFilter.h"
#include "itkCyclicShiftImageFilter.h"
#include "itkExtractImageFilter.h"
#include "itkFFTPadImageFilter.h"
#include "itkImageBase.h"
#include "itkMultiplyImageFilter.h"
#include "itkNormalizeToConstantImageFilter.h"
#include "itkMath.h"
#include "itkRegionOfInterestImageFilter.h"

namespace itk
{

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::FFTConvolutionImageFilter()
{
  m_SizeGreatestPrimeFactor = FFTFilterType::New()->GetSizeGreatestPrimeFactor();
}

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
void
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::GenerateInputRequestedRegion()
{
  // Pad the input image with the radius of the kernel.
  if (this->GetInput())
  {
    InputRegionType inputRegion = this->GetOutput()->GetRequestedRegion();

    // Pad the output request region by the kernel radius.
    // Note that FFT padding will always be generated by a boundary condition
    // as an implementation detail, while pixels for kernel radius padding may be taken
    // from the original image if they lies inside the image bounds.
    inputRegion.PadByRadius(this->GetKernelRadius());

    // Crop the output requested region to fit within the largest
    // possible region.
    InputImageType * inputPtr = itkDynamicCastInDebugMode<InputImageType *>(this->GetPrimaryInput());
    bool             wasPartiallyInside = inputRegion.Crop(inputPtr->GetLargestPossibleRegion());
    if (!wasPartiallyInside)
    {
      itkExceptionMacro("Requested region is outside the largest possible region.");
    }

    // Input is an image, cast away the constness so we can set
    // the requested region.
    inputPtr->SetRequestedRegion(inputRegion);
  }

  if (this->GetKernelImage())
  {
    // Input kernel is an image, cast away the constness so we can set
    // the requested region.
    typename KernelImageType::Pointer kernelPtr = const_cast<KernelImageType *>(this->GetKernelImage());
    kernelPtr->SetRequestedRegionToLargestPossibleRegion();
  }
}

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
void
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::GenerateData()
{
  // Create a process accumulator for tracking the progress of this minipipeline
  auto progress = ProgressAccumulator::New();
  progress->SetMiniPipelineFilter(this);

  auto localInput = InputImageType::New();
  localInput->Graft(this->GetInput());

  const KernelImageType * kernelImage = this->GetKernelImage();

  InternalComplexImagePointerType input = nullptr;
  InternalComplexImagePointerType kernel = nullptr;
  this->PrepareInputs(localInput, kernelImage, input, kernel, progress, 0.7f);

  // Convolve
  using MultiplyFilterType =
    MultiplyImageFilter<InternalComplexImageType, InternalComplexImageType, InternalComplexImageType>;
  auto multiplyFilter = MultiplyFilterType::New();
  multiplyFilter->SetInput1(input);
  multiplyFilter->SetInput2(kernel);
  multiplyFilter->ReleaseDataFlagOn();
  progress->RegisterInternalFilter(multiplyFilter, 0.1);
  multiplyFilter->Update();

  // Free up the memory for the prepared inputs
  input = nullptr;
  kernel = nullptr;

  this->ProduceOutput(multiplyFilter->GetOutput(), progress, 0.2);
}

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
void
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::PrepareInputs(
  const InputImageType *            input,
  const KernelImageType *           kernel,
  InternalComplexImagePointerType & preparedInput,
  InternalComplexImagePointerType & preparedKernel,
  ProgressAccumulator *             progress,
  float                             progressWeight)
{
  this->PrepareInput(input, preparedInput, progress, 0.5f * progressWeight);
  this->PrepareKernel(kernel, preparedKernel, progress, 0.5f * progressWeight);
}

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
void
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::PrepareInput(
  const InputImageType *            input,
  InternalComplexImagePointerType & preparedInput,
  ProgressAccumulator *             progress,
  float                             progressWeight)
{
  InternalImagePointerType paddedInput;
  this->PadInput(input, paddedInput, progress, 0.3f * progressWeight);
  this->TransformPaddedInput(paddedInput, preparedInput, progress, 0.7f * progressWeight);
}

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
void
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::PadInput(
  const InputImageType *     input,
  InternalImagePointerType & paddedInput,
  ProgressAccumulator *      progress,
  float                      progressWeight)
{
  float remainingProgress = 1.0f;

  InputRegionType  inputLargestRegion = input->GetLargestPossibleRegion();
  InputSizeType    inputLargestSize = inputLargestRegion.GetSize();
  InputIndexType   inputLargestIndex = inputLargestRegion.GetIndex();
  InputRegionType  inputRequestedRegion = input->GetRequestedRegion();
  InputSizeType    inputRequestedSize = inputRequestedRegion.GetSize();
  InputIndexType   inputRequestedIndex = inputRequestedRegion.GetIndex();
  OutputRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
  OutputSizeType   outputRequestedSize = outputRequestedRegion.GetSize();
  OutputIndexType  outputRequestedIndex = outputRequestedRegion.GetIndex();

  // Pad the input image such that the requested region, expanded by
  // twice the kernel radius, lies entirely within the buffered region.
  // If the requested region is at least the kernel radius away from
  // all boundaries of the buffered input region then padding is skipped.

  using InputPadFilterType = PadImageFilter<InputImageType, InputImageType>;
  using PadSizeType = typename InputPadFilterType::SizeType;

  PadSizeType    lowerPad, upperPad;
  bool           needsKernelPadding = false;
  KernelSizeType kernelRadius = this->GetKernelRadius();
  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
  {
    // Use signed int for arithmetic
    int largestUpperCorner = inputLargestIndex[dim] + inputLargestSize[dim];
    int requestedUpperCorner = inputRequestedIndex[dim] + inputRequestedSize[dim];

    int largestLowerCorner = inputLargestIndex[dim];
    int requestedLowerCorner = inputRequestedIndex[dim];

    // Pad for difference between lower corner of largest vs requested region
    int lower = static_cast<int>(kernelRadius[dim]) - (requestedLowerCorner - largestLowerCorner);
    // Pad for difference between upper corner of largest vs requested region
    int upper = static_cast<int>(kernelRadius[dim]) - (largestUpperCorner - requestedUpperCorner);
    lowerPad[dim] = (lower > 0 ? lower : 0);
    upperPad[dim] = (upper > 0 ? upper : 0);

    needsKernelPadding = needsKernelPadding || lower > 0 || upper > 0;
  }

  const InputImageType * kernelPaddedInput;
  if (needsKernelPadding)
  {
    auto inputPadder = InputPadFilterType::New();
    inputPadder->SetBoundaryCondition(this->GetBoundaryCondition());
    inputPadder->SetPadLowerBound(lowerPad);
    inputPadder->SetPadUpperBound(upperPad);
    inputPadder->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
    inputPadder->SetInput(input);
    inputPadder->ReleaseDataFlagOn();
    progress->RegisterInternalFilter(inputPadder, 0.2f * progressWeight);
    remainingProgress -= 0.2f;
    inputPadder->Update();
    kernelPaddedInput = inputPadder->GetOutput();
  }
  else
  {
    kernelPaddedInput = input;
  }

  // Crop to region of interest to minimize FFT requested region size.
  // This blocks forward FFT from requesting the entire input image
  // when only a subregion (output requested region + kernel radius padding)
  // is required for convolution to succeed.
  // TODO: Improve readability by keeping the output index as-is rather than resetting to 0, if the filter allows it

  const InputImageType * regionOfInterestImage;
  InputRegionType        regionOfInterest;
  InputSizeType          regionOfInterestSize;
  InputIndexType         regionOfInterestIndex;

  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
  {
    regionOfInterestSize[dim] = outputRequestedSize[dim] + 2 * kernelRadius[dim];
    regionOfInterestIndex[dim] = outputRequestedIndex[dim] - kernelRadius[dim];
  }
  regionOfInterest = InputRegionType(regionOfInterestIndex, regionOfInterestSize);

  if (outputRequestedRegion != inputLargestRegion)
  {

    using CropFilterType = RegionOfInterestImageFilter<InputImageType, InputImageType>;
    auto cropFilter = CropFilterType::New();
    cropFilter->SetInput(kernelPaddedInput);
    cropFilter->SetRegionOfInterest(regionOfInterest);
    cropFilter->ReleaseDataFlagOn();
    progress->RegisterInternalFilter(cropFilter, 0.1f * progressWeight);
    remainingProgress -= 0.1f;

    // Shift region of interest back to original index
    // so that filter output is physically correct
    using InfoFilterType = ChangeInformationImageFilter<InputImageType>;
    using InfoOffsetValueType = typename InfoFilterType::OutputImageOffsetValueType;
    auto inputInfoFilter = InfoFilterType::New();
    inputInfoFilter->SetInput(cropFilter->GetOutput());
    inputInfoFilter->ChangeRegionOn();
    InfoOffsetValueType offset[ImageDimension];
    for (unsigned int dim = 0; dim < ImageDimension; ++dim)
    {
      offset[dim] = regionOfInterestIndex[dim];
    }
    inputInfoFilter->SetOutputOffset(offset);
    inputInfoFilter->ReleaseDataFlagOn();
    inputInfoFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
    progress->RegisterInternalFilter(inputInfoFilter, 0.001f * progressWeight);
    remainingProgress -= 0.001f;
    inputInfoFilter->Update();
    regionOfInterestImage = inputInfoFilter->GetOutput();
  }
  else
  {
    regionOfInterestImage = kernelPaddedInput;
  }

  // Pad for FFT so that each image side is factorable by at most
  // the specified greatest prime factor.
  // Note that FFT padding is always taken from the boundary condition
  // to avoid introducing extra information to FFT convolution vs spatial convolution
  using FFTPadFilterType = itk::FFTPadImageFilter<InputImageType, InputImageType>;
  auto fftPadder = FFTPadFilterType::New();
  fftPadder->SetInput(regionOfInterestImage);
  fftPadder->SetSizeGreatestPrimeFactor(m_SizeGreatestPrimeFactor);
  fftPadder->SetBoundaryCondition(this->GetBoundaryCondition());
  fftPadder->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
  fftPadder->ReleaseDataFlagOn();
  progress->RegisterInternalFilter(fftPadder, 0.199f * progressWeight);
  remainingProgress -= 0.199f;
  fftPadder->Update();

  // Save FFT padding size to assist in output cropping
  InputSizeType fftPaddedSize = fftPadder->GetOutput()->GetLargestPossibleRegion().GetSize();
  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
  {
    m_FFTPadSize[dim] = fftPaddedSize[dim] - regionOfInterestSize[dim];
  }

  // We could avoid a separate cast here by setting the output type of
  // the padder to the InternalImageType, but doing so complicates the
  // definition of the boundary condition passed into this class and
  // requires the InternalImageType to be exposed publicly.
  // TODO: Casting has been observed to take the bulk of image prep time, up to 75%.
  //       Optimize to avoid unnecessary casting and copying of memory.
  using InputCastFilterType = CastImageFilter<InputImageType, InternalImageType>;
  auto inputCaster = InputCastFilterType::New();
  inputCaster->InPlaceOn();
  inputCaster->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
  inputCaster->SetInput(fftPadder->GetOutput());
  inputCaster->ReleaseDataFlagOn();
  progress->RegisterInternalFilter(inputCaster, remainingProgress * progressWeight);
  inputCaster->Update();

  // Cache padded input size to use in kernel preparation and output cropping
  m_PaddedInputRegion = inputCaster->GetOutput()->GetLargestPossibleRegion();
  paddedInput = inputCaster->GetOutput();
}

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
void
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::TransformPaddedInput(
  const InternalImageType *         paddedInput,
  InternalComplexImagePointerType & transformedInput,
  ProgressAccumulator *             progress,
  float                             progressWeight)
{
  // Take the Fourier transform of the padded image.
  auto imageFFTFilter = FFTFilterType::New();
  imageFFTFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
  imageFFTFilter->SetInput(paddedInput);
  imageFFTFilter->ReleaseDataFlagOn();
  progress->RegisterInternalFilter(imageFFTFilter, progressWeight);
  imageFFTFilter->Update();

  transformedInput = imageFFTFilter->GetOutput();
  transformedInput->DisconnectPipeline();

  imageFFTFilter->SetInput(nullptr);
  imageFFTFilter = nullptr;
}

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
void
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::PrepareKernel(
  const KernelImageType *           kernel,
  InternalComplexImagePointerType & preparedKernel,
  ProgressAccumulator *             progress,
  float                             progressWeight)
{
  KernelRegionType kernelRegion = kernel->GetLargestPossibleRegion();
  KernelSizeType   kernelSize = kernelRegion.GetSize();

  InputSizeType                      inputPadSize = m_PaddedInputRegion.GetSize();
  typename KernelImageType::SizeType kernelUpperBound;
  for (unsigned int i = 0; i < ImageDimension; ++i)
  {
    kernelUpperBound[i] = inputPadSize[i] - kernelSize[i];
  }

  InternalImagePointerType paddedKernelImage = nullptr;

  float paddingWeight = 0.2f;
  if (this->GetNormalize())
  {
    using NormalizeFilterType = NormalizeToConstantImageFilter<KernelImageType, InternalImageType>;
    auto normalizeFilter = NormalizeFilterType::New();
    normalizeFilter->SetConstant(NumericTraits<TInternalPrecision>::OneValue());
    normalizeFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
    normalizeFilter->SetInput(kernel);
    normalizeFilter->ReleaseDataFlagOn();
    progress->RegisterInternalFilter(normalizeFilter, 0.2f * paddingWeight * progressWeight);

    // Pad the kernel image with zeros.
    using KernelPadType = ConstantPadImageFilter<InternalImageType, InternalImageType>;
    using KernelPadPointer = typename KernelPadType::Pointer;
    KernelPadPointer kernelPadder = KernelPadType::New();
    kernelPadder->SetConstant(TInternalPrecision{});
    kernelPadder->SetPadUpperBound(kernelUpperBound);
    kernelPadder->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
    kernelPadder->SetInput(normalizeFilter->GetOutput());
    kernelPadder->ReleaseDataFlagOn();
    progress->RegisterInternalFilter(kernelPadder, 0.8f * paddingWeight * progressWeight);
    kernelPadder->Update();
    paddedKernelImage = kernelPadder->GetOutput();
  }
  else
  {
    // Pad the kernel image with zeros.
    using KernelPadType = ConstantPadImageFilter<KernelImageType, InternalImageType>;
    using KernelPadPointer = typename KernelPadType::Pointer;
    KernelPadPointer kernelPadder = KernelPadType::New();
    kernelPadder->SetConstant(TInternalPrecision{});
    kernelPadder->SetPadUpperBound(kernelUpperBound);
    kernelPadder->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
    kernelPadder->SetInput(kernel);
    kernelPadder->ReleaseDataFlagOn();
    progress->RegisterInternalFilter(kernelPadder, paddingWeight * progressWeight);
    paddedKernelImage = kernelPadder->GetOutput();
  }

  // Shift the padded kernel image.
  using KernelShiftFilterType = CyclicShiftImageFilter<InternalImageType, InternalImageType>;
  auto                                       kernelShifter = KernelShiftFilterType::New();
  typename KernelShiftFilterType::OffsetType kernelShift;
  for (unsigned int i = 0; i < ImageDimension; ++i)
  {
    kernelShift[i] = -(static_cast<typename KernelShiftFilterType::OffsetType::OffsetValueType>(kernelSize[i] / 2));
  }
  kernelShifter->SetShift(kernelShift);
  kernelShifter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
  kernelShifter->SetInput(paddedKernelImage);
  kernelShifter->ReleaseDataFlagOn();
  progress->RegisterInternalFilter(kernelShifter, 0.1f * progressWeight);

  // Compute the kernel complex image
  auto kernelFFTFilter = FFTFilterType::New();
  kernelFFTFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
  kernelFFTFilter->SetInput(kernelShifter->GetOutput());
  progress->RegisterInternalFilter(kernelFFTFilter, 0.699f * progressWeight);
  kernelFFTFilter->Update();

  // Shift the kernel complex image in space so that it coincides with the
  // input complex image
  using InfoFilterType = ChangeInformationImageFilter<InternalComplexImageType>;
  auto kernelInfoFilter = InfoFilterType::New();
  kernelInfoFilter->ChangeRegionOn();

  using InfoOffsetValueType = typename InfoFilterType::OutputImageOffsetValueType;
  const InputIndexType &  inputIndex = m_PaddedInputRegion.GetIndex();
  const KernelIndexType & kernelIndex = kernel->GetLargestPossibleRegion().GetIndex();
  InfoOffsetValueType     kernelOffset[ImageDimension];
  for (unsigned int i = 0; i < ImageDimension; ++i)
  {
    kernelOffset[i] = static_cast<InfoOffsetValueType>(inputIndex[i] - kernelIndex[i]);
  }
  kernelInfoFilter->SetOutputOffset(kernelOffset);
  kernelInfoFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
  kernelInfoFilter->SetInput(kernelFFTFilter->GetOutput());
  progress->RegisterInternalFilter(kernelInfoFilter, 0.001f * progressWeight);
  kernelInfoFilter->Update();

  preparedKernel = kernelInfoFilter->GetOutput();
}

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
void
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::ProduceOutput(
  InternalComplexImageType * paddedOutput,
  ProgressAccumulator *      progress,
  float                      progressWeight)
{
  // Transform the convolution result back into physical space
  auto ifftFilter = IFFTFilterType::New();
  ifftFilter->SetActualXDimensionIsOdd(this->GetXDimensionIsOdd());
  ifftFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
  ifftFilter->SetInput(paddedOutput);
  ifftFilter->ReleaseDataFlagOn();
  progress->RegisterInternalFilter(ifftFilter, 0.6f * progressWeight);
  ifftFilter->Update();

  this->CropOutput(ifftFilter->GetOutput(), progress, 0.4f * progressWeight);
}

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
void
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::CropOutput(
  InternalImageType *   paddedOutput,
  ProgressAccumulator * progress,
  float                 progressWeight)
{
  InternalIndexType paddedOutputIndex = paddedOutput->GetLargestPossibleRegion().GetIndex();

  // Allocate the output
  this->AllocateOutputs();

  // Now crop the output to the desired size.
  using ExtractFilterType = ExtractImageFilter<InternalImageType, OutputImageType>;

  auto extractFilter = ExtractFilterType::New();
  extractFilter->InPlaceOn();
  extractFilter->GraftOutput(this->GetOutput());

  // Set up the crop sizes.
  InternalIndexType extractionIndex;
  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
  {
    extractionIndex[dim] = paddedOutputIndex[dim] + m_FFTPadSize[dim] / 2 + GetKernelRadius()[dim];
  }
  auto               requestedRegionSize = this->GetOutput()->GetRequestedRegion().GetSize();
  InternalRegionType extractionRegion(extractionIndex, requestedRegionSize);
  extractFilter->SetExtractionRegion(extractionRegion);

  // Graft the minipipeline output to this filter.
  extractFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
  extractFilter->SetInput(paddedOutput);
  extractFilter->GetOutput()->SetRequestedRegion(this->GetOutput()->GetRequestedRegion());
  progress->RegisterInternalFilter(extractFilter, progressWeight);
  extractFilter->Update();

  OutputImageType * extractedImage = extractFilter->GetOutput();
  OutputImageType * output = this->GetOutput();

  // Only manually copy the buffer via the pixel container.
  // The output meta-data of the extract filter is not correct and
  // different that the GenerateOutputInformation method. So just copy
  // the buffer.
  output->SetBufferedRegion(extractedImage->GetBufferedRegion());
  output->SetPixelContainer(extractedImage->GetPixelContainer());
}

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
auto
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::GetKernelRadius() const
  -> KernelSizeType
{
  KernelSizeType kernelSize = this->GetKernelImage()->GetLargestPossibleRegion().GetSize();
  KernelSizeType kernelRadius;
  for (unsigned int i = 0; i < ImageDimension; ++i)
  {
    kernelRadius[i] = kernelSize[i] / 2;
  }
  return kernelRadius;
}

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
bool
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::GetXDimensionIsOdd() const
{
  return (m_PaddedInputRegion.GetSize()[0] % 2 != 0);
}

template <typename TInputImage, typename TKernelImage, typename TOutputImage, typename TInternalPrecision>
void
FFTConvolutionImageFilter<TInputImage, TKernelImage, TOutputImage, TInternalPrecision>::PrintSelf(std::ostream & os,
                                                                                                  Indent indent) const
{
  Superclass::PrintSelf(os, indent);
  os << indent << "SizeGreatestPrimeFactor: " << m_SizeGreatestPrimeFactor << std::endl;
}

} // namespace itk
#endif // itkFFTConvolutionImageFilter_hxx