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

#include "itkPoint.h"
#include "itkMath.h"

namespace itk
{
template <unsigned int VDimension>
auto
BresenhamLine<VDimension>::BuildLine(LType Direction, IdentifierType length) -> OffsetArray
{
  // copied from the line iterator
  /** Variables that drive the Bresenham-Algorithm */
  // The dimension with the largest difference between start and end
  unsigned int m_MainDirection;

  // Accumulated error for the other dimensions
  IndexType m_AccumulateError;

  // Increment for the error for each step. Two times the difference between
  // start and end
  IndexType m_IncrementError;

  // If enough is accumulated for a dimension, the index has to be
  // incremented. Will be the number of pixels in the line
  IndexType m_MaximalError;

  // Direction of increment. -1 or 1
  IndexType m_OverflowIncrement;

  // After an overflow, the accumulated error is reduced again. Will be
  // two times the number of pixels in the line
  IndexType m_ReduceErrorAfterIncrement;

  OffsetArray result(length);

  IndexType m_CurrentImageIndex, LastIndex;

  Direction.Normalize();
  // we are going to start at 0
  m_CurrentImageIndex.Fill(0);
  constexpr IndexType StartIndex = { { 0 } };
  for (unsigned int i = 0; i < VDimension; ++i)
  {
    LastIndex[i] = (IndexValueType)(length * Direction[i]);
  }
  // Find the dominant direction
  IndexValueType maxDistance = 0;
  unsigned int   maxDistanceDimension = 0;
  for (unsigned int i = 0; i < VDimension; ++i)
  {
    auto distance = static_cast<long>(itk::Math::abs(LastIndex[i]));
    if (distance > maxDistance)
    {
      maxDistance = distance;
      maxDistanceDimension = i;
    }
    m_IncrementError[i] = 2 * distance;
    m_OverflowIncrement[i] = (LastIndex[i] < 0 ? -1 : 1);
  }
  m_MainDirection = maxDistanceDimension;
  m_MaximalError.Fill(maxDistance);
  m_ReduceErrorAfterIncrement.Fill(2 * maxDistance);
  m_AccumulateError.Fill(0);
  unsigned int steps = 1;
  result[0] = m_CurrentImageIndex - StartIndex;
  while (steps < length)
  {
    // This part is from ++ in LineConstIterator
    // We need to modify m_AccumulateError, m_CurrentImageIndex, m_IsAtEnd
    for (unsigned int i = 0; i < VDimension; ++i)
    {
      if (i == m_MainDirection)
      {
        m_CurrentImageIndex[i] += m_OverflowIncrement[i];
      }
      else
      {
        m_AccumulateError[i] += m_IncrementError[i];
        if (m_AccumulateError[i] >= m_MaximalError[i])
        {
          m_CurrentImageIndex[i] += m_OverflowIncrement[i];
          m_AccumulateError[i] -= m_ReduceErrorAfterIncrement[i];
        }
      }
    }

    result[steps] = m_CurrentImageIndex - StartIndex; // produce an offset

    ++steps;
  }
  return (result);
}

template <unsigned int VDimension>
auto
BresenhamLine<VDimension>::BuildLine(IndexType p0, IndexType p1) -> IndexArray
{
  itk::Point<float, VDimension> point0;
  itk::Point<float, VDimension> point1;
  IdentifierType                maxDistance = 0;
  for (unsigned int i = 0; i < VDimension; ++i)
  {
    point0[i] = p0[i];
    point1[i] = p1[i];
    IdentifierType distance = itk::Math::abs(p0[i] - p1[i]) + 1;
    if (distance > maxDistance)
    {
      maxDistance = distance;
    }
  }

  OffsetArray offsets = this->BuildLine(point1 - point0, maxDistance + 1);

  IndexArray indices;
  indices.reserve(offsets.size()); // we might not have to use the last one
  for (unsigned int i = 0; i < offsets.size(); ++i)
  {
    indices.push_back(p0 + offsets[i]);
    if (indices.back() == p1)
    {
      break;
    }
  }
  return indices;
}

} // namespace itk

#endif