/*=========================================================================
 *
 *  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 itkGaborImageSource_hxx
#define itkGaborImageSource_hxx

#include "itkGaborKernelFunction.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkProgressReporter.h"
#include "itkObjectFactory.h"

namespace itk
{

template <typename TOutputImage>
GaborImageSource<TOutputImage>::GaborImageSource()

{
  // Gabor parameters, defined so that the Gaussian
  // is centered in the default image
  this->m_Mean.Fill(32.0);
  this->m_Sigma.Fill(16.0);
}

template <typename TOutputImage>
void
GaborImageSource<TOutputImage>::GenerateData()
{
  OutputImageType * outputPtr = this->GetOutput();

  // Allocate the output buffer
  outputPtr->SetBufferedRegion(outputPtr->GetRequestedRegion());
  outputPtr->Allocate();

  // Create and initialize a new Gaussian function
  using KernelFunctionType = GaborKernelFunction<double>;
  auto gabor = KernelFunctionType::New();

  gabor->SetSigma(this->m_Sigma[0]);
  gabor->SetFrequency(this->m_Frequency);
  gabor->SetPhaseOffset(this->m_PhaseOffset);
  gabor->SetCalculateImaginaryPart(this->m_CalculateImaginaryPart);

  ProgressReporter progress(this, 0, outputPtr->GetRequestedRegion().GetNumberOfPixels());

  // Walk the output image, evaluating the spatial function at each pixel
  for (ImageRegionIteratorWithIndex<OutputImageType> outIt(outputPtr, outputPtr->GetRequestedRegion());
       !outIt.IsAtEnd();
       ++outIt)
  {
    const typename OutputImageType::IndexType index = outIt.GetIndex();
    // The position at which the function is evaluated
    typename OutputImageType::PointType evalPoint;
    outputPtr->TransformIndexToPhysicalPoint(index, evalPoint);
    double sum = 0.0;
    for (unsigned int i = 1; i < ImageDimension; ++i)
    {
      sum += itk::Math::sqr((evalPoint[i] - this->m_Mean[i]) / this->m_Sigma[i]);
    }
    const double value = std::exp(-0.5 * sum) * gabor->Evaluate(evalPoint[0] - this->m_Mean[0]);

    // Set the pixel value to the function value
    outIt.Set(static_cast<PixelType>(value));
    progress.CompletedPixel();
  }
}

template <typename TOutputImage>
void
GaborImageSource<TOutputImage>::PrintSelf(std::ostream & os, Indent indent) const
{
  Superclass::PrintSelf(os, indent);

  if (this->GetCalculateImaginaryPart())
  {
    os << indent << "Calculate complex part: true " << std::endl;
  }
  else
  {
    os << indent << "Calculate complex part: false " << std::endl;
  }
  os << indent << "Frequency: " << this->GetFrequency() << std::endl;
  os << indent << "Phase offset: " << m_PhaseOffset << std::endl;
  os << indent << "Sigma: " << this->GetSigma() << std::endl;
  os << indent << "Mean: " << this->GetMean() << std::endl;
}
} // end namespace itk

#endif