/*=========================================================================
 *
 *  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
 *
 *         http://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 itkImageKernelOperator_hxx
#define itkImageKernelOperator_hxx

#include "itkImageKernelOperator.h"

#include "itkImageBufferRange.h"
#include "itkImageRegionConstIterator.h"

/*
 *
 * This code was contributed in the Insight Journal paper:
 *
 * "Image Kernel Convolution"
 * by Tustison N., Gee J.
 * https://hdl.handle.net/1926/1323
 * http://www.insight-journal.org/browse/publication/208
 *
 */

namespace itk
{

template <typename TPixel, unsigned int VDimension, typename TAllocator>
void
ImageKernelOperator<TPixel, VDimension, TAllocator>::SetImageKernel(const ImageType * kernel)
{
  m_ImageKernel = kernel;
}

template <typename TPixel, unsigned int VDimension, typename TAllocator>
const typename ImageKernelOperator<TPixel, VDimension, TAllocator>::ImageType *
ImageKernelOperator<TPixel, VDimension, TAllocator>::GetImageKernel() const
{
  return m_ImageKernel;
}

template <typename TPixel, unsigned int VDimension, typename TAllocator>
typename ImageKernelOperator<TPixel, VDimension, TAllocator>::CoefficientVector
ImageKernelOperator<TPixel, VDimension, TAllocator>::GenerateCoefficients()
{
  // Check that the input image is fully buffered.
  if (m_ImageKernel->GetBufferedRegion() != m_ImageKernel->GetLargestPossibleRegion())
  {
    itkExceptionMacro(<< "ImageKernel is not fully buffered. " << std::endl
                      << "Buffered region: " << m_ImageKernel->GetBufferedRegion() << std::endl
                      << "Largest possible region: " << m_ImageKernel->GetLargestPossibleRegion() << std::endl
                      << "You should call UpdateLargestPossibleRegion() on "
                      << "the filter whose output is passed to "
                      << "SetImageKernel().");
  }

  // Check that the size of the kernel is odd in all dimensions.
  for (unsigned int i = 0; i < VDimension; i++)
  {
    if (m_ImageKernel->GetLargestPossibleRegion().GetSize()[i] % 2 == 0)
    {
      itkExceptionMacro(<< "ImageKernelOperator requires an input image "
                        << "whose size is odd in all dimensions. The provided "
                        << "image has size " << m_ImageKernel->GetLargestPossibleRegion().GetSize());
    }
  }

  const auto imageBufferRange = Experimental::MakeImageBufferRange(m_ImageKernel.GetPointer());

  return CoefficientVector(imageBufferRange.cbegin(), imageBufferRange.cend());
}

template <typename TPixel, unsigned int VDimension, typename TAllocator>
void
ImageKernelOperator<TPixel, VDimension, TAllocator>::Fill(const CoefficientVector & coeff)
{
  // Initialize all coefficients to zero
  this->InitializeToZero();

  std::slice *                               temp_slice;
  typename CoefficientVector::const_iterator it;

  temp_slice = new std::slice(0, coeff.size(), 1);
  it = coeff.begin();

  typename Superclass::SliceIteratorType data(this, *temp_slice);
  delete temp_slice;

  // Copy the coefficients into the neighborhood, truncating them if there
  // are too many.
  for (data = data.Begin(); data < data.End(); ++data, ++it)
  {
    *data = static_cast<TPixel>(*it);
  }
}
} // end namespace itk

#endif
