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


namespace itk
{
template <typename TInputImage, typename TOutputImage>
auto
ZeroFluxNeumannBoundaryCondition<TInputImage, TOutputImage>::operator()(const OffsetType &       point_index,
                                                                        const OffsetType &       boundary_offset,
                                                                        const NeighborhoodType * data) const
  -> OutputPixelType
{
  int linear_index = 0;

  // Return the value of the pixel at the closest boundary point.
  for (unsigned int i = 0; i < ImageDimension; ++i)
  {
    linear_index += (point_index[i] + boundary_offset[i]) * data->GetStride(i);
  }

  // The reinterpret_cast is necessary, cause we will have a warning if we
  // do not do this. (In fact this function exists for legacy
  // reasons. The overloaded function below should be (and is) used instead).
  // See any of the neighborhood iterators.
  //
  // (data->operator[](linear_index)) is guaranteed to be a pointer to
  // TImage::PixelType except for VectorImage, in which case, it will be a
  // pointer to TImage::InternalPixelType.
  //
  // A typical neighborhood iterator working on an image will use the boundary
  // condition in the following manner:
  //
  // \code
  // // Initialize the functor typically in the constructor.
  // m_NeighborhoodAccessorFunctor = image->GetNeighborhoodAccessor();
  // m_NeighborhoodAccessorFunctor->SetBegin( image->GetBufferPointer() );
  //
  // m_NeighborhoodAccessorFunctor.BoundaryCondition(
  //    point_index, boundary_offset, data, m_ChosenBoundaryCondition );
  // \endcode
  //
  return static_cast<OutputPixelType>(*(reinterpret_cast<PixelType *>((data->operator[](linear_index)))));
}

template <typename TInputImage, typename TOutputImage>
auto
ZeroFluxNeumannBoundaryCondition<TInputImage, TOutputImage>::operator()(
  const OffsetType &                      point_index,
  const OffsetType &                      boundary_offset,
  const NeighborhoodType *                data,
  const NeighborhoodAccessorFunctorType & neighborhoodAccessorFunctor) const -> OutputPixelType
{
  int linear_index = 0;

  // Return the value of the pixel at the closest boundary point.
  for (unsigned int i = 0; i < ImageDimension; ++i)
  {
    linear_index += (point_index[i] + boundary_offset[i]) * data->GetStride(i);
  }

  return static_cast<OutputPixelType>(neighborhoodAccessorFunctor.Get(data->operator[](linear_index)));
}


template <typename TInputImage, typename TOutputImage>
auto
ZeroFluxNeumannBoundaryCondition<TInputImage, TOutputImage>::GetInputRequestedRegion(
  const RegionType & inputLargestPossibleRegion,
  const RegionType & outputRequestedRegion) const -> RegionType
{
  IndexType inputIndex = inputLargestPossibleRegion.GetIndex();
  SizeType  inputSize = inputLargestPossibleRegion.GetSize();

  IndexType outputIndex = outputRequestedRegion.GetIndex();
  SizeType  outputSize = outputRequestedRegion.GetSize();

  IndexType requestIndex;
  SizeType  requestSize;

  for (unsigned int i = 0; i < ImageDimension; ++i)
  {
    // Check if the output region is entirely below the low index of
    // the image region.
    if (outputIndex[i] + static_cast<OffsetValueType>(outputSize[i]) <= inputIndex[i])
    {
      // Include an image layer one pixel thick closest to the outputRequestedRegion
      requestIndex[i] = inputIndex[i];
      requestSize[i] = 1;
    }
    // Check if the output is entirely above the high index of the
    // image region.
    else if (outputIndex[i] >= inputIndex[i] + static_cast<OffsetValueType>(inputSize[i]))
    {
      // Include an image layer one pixel thick closest to the outputRequestedRegion
      requestIndex[i] = inputIndex[i] + static_cast<OffsetValueType>(inputSize[i]) - 1;
      requestSize[i] = 1;
    }
    // The output region intersects the image region.
    else
    {
      requestIndex[i] = inputIndex[i];
      requestSize[i] = inputSize[i];

      // First check the start index
      if (requestIndex[i] < outputIndex[i])
      {
        // How much do we need to adjust
        OffsetValueType crop = outputIndex[i] - requestIndex[i];

        // Adjust the start index and the size of the current region
        requestIndex[i] += crop;
        requestSize[i] -= static_cast<SizeValueType>(crop);
      }
      // Now check the final size
      if (requestIndex[i] + static_cast<OffsetValueType>(requestSize[i]) >
          outputIndex[i] + static_cast<OffsetValueType>(outputSize[i]))
      {
        // How much do we need to adjust
        OffsetValueType crop = requestIndex[i] + static_cast<OffsetValueType>(requestSize[i]) - outputIndex[i] -
                               static_cast<OffsetValueType>(outputSize[i]);

        // Adjust the size
        requestSize[i] -= static_cast<SizeValueType>(crop);
      }
    }
  }

  RegionType inputRequestedRegion(requestIndex, requestSize);

  return inputRequestedRegion;
}


template <typename TInputImage, typename TOutputImage>
auto
ZeroFluxNeumannBoundaryCondition<TInputImage, TOutputImage>::GetPixel(const IndexType &   index,
                                                                      const TInputImage * image) const
  -> OutputPixelType
{
  RegionType imageRegion = image->GetLargestPossibleRegion();
  IndexType  imageIndex = imageRegion.GetIndex();
  SizeType   imageSize = imageRegion.GetSize();

  IndexType lookupIndex;

  for (unsigned int i = 0; i < ImageDimension; ++i)
  {
    IndexValueType lowerIndex = imageIndex[i];
    IndexValueType upperIndex = imageIndex[i] + static_cast<IndexValueType>(imageSize[i]) - 1;
    if (index[i] < lowerIndex)
    {
      lookupIndex[i] = lowerIndex;
    }
    else if (index[i] > upperIndex)
    {
      lookupIndex[i] = upperIndex;
    }
    else // in bounds
    {
      lookupIndex[i] = index[i];
    }
  }

  return static_cast<OutputPixelType>(image->GetPixel(lookupIndex));
}

} // end namespace itk

#endif