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

#include "itkImageRegionConstIteratorWithIndex.h"
#include "itkImageRegionConstIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include <list>

namespace itk
{
/**
 * \class SharedMorphUtilities
 * \brief functionality in common for anchor and VanHerkGilWerman openings/closings and
 * erosions/dilation
 *
 */

template <typename TRegion, typename TLine>
bool
NeedToDoFace(const TRegion AllImage, const TRegion face, const TLine line)
{
  // can't use the continuous IsInside (even if I could get it to
  // work) because on the edge doesn't count as inside for this test

  // If the component of the vector orthogonal to the face doesn't go
  // inside the image then we can ignore the face

  // find the small dimension of the face - should only be one
  typename TRegion::IndexType ISt = AllImage.GetIndex();

  typename TRegion::SizeType  FSz = face.GetSize();
  typename TRegion::IndexType FSt = face.GetIndex();

  unsigned int smallDim = 0;
  for (unsigned int i = 0; i < AllImage.GetImageDimension(); ++i)
  {
    if (FSz[i] == 1)
    {
      smallDim = i;
      break;
    }
  }
  IndexValueType startI = ISt[smallDim];
  IndexValueType facePos = FSt[smallDim] + FSz[smallDim] - 1;
  if (facePos == startI)
  {
    // at the start of dimension - vector must be positive
    if (line[smallDim] > 0.000001)
    {
      return true;
    }
    // some small angle that we consider to be zero - should be more rigorous
  }
  else
  {
    // at the end of dimension - vector must be positive
    if (line[smallDim] < -0.000001)
    {
      return true;
    }
  }
  return (false);
}

template <typename TImage, typename TBres, typename TLine>
int
ComputeStartEnd(const typename TImage::IndexType  StartIndex,
                const TLine                       line,
                const float                       tol,
                const typename TBres::OffsetArray LineOffsets,
                const typename TImage::RegionType AllImage,
                unsigned int &                    start,
                unsigned int &                    end)
{
  // compute intersection between ray and box
  typename TImage::IndexType ImStart = AllImage.GetIndex();
  typename TImage::SizeType  ImSize = AllImage.GetSize();
  float                      Tfar = NumericTraits<float>::max();
  float                      Tnear = NumericTraits<float>::NonpositiveMin();
  float                      domdir = NumericTraits<float>::NonpositiveMin();
  int                        sPos, ePos;
  unsigned int               perpdir = 0;
  for (unsigned int i = 0; i < TImage::RegionType::ImageDimension; ++i)
  {
    const auto abs_line_elmt_tmp = itk::Math::abs(line[i]);
    if (abs_line_elmt_tmp > domdir)
    {
      domdir = abs_line_elmt_tmp;
      perpdir = i;
    }
    if (abs_line_elmt_tmp > tol)
    {
      const float P1{ static_cast<float>(ImStart[i] - StartIndex[i]) };
      const float P2{ P1 + static_cast<float>(ImSize[i]) - 1.0F };
      float       T1 = P1 / line[i];
      float       T2 = P2 / line[i];

      if (T1 > T2)
      {
        // T1 is meant to be the near face
        std::swap(T1, T2);
      }
      // want the farthest Tnear and nearest TFar
      if (T1 > Tnear)
      {
        Tnear = T1;
      }
      if (T2 < Tfar)
      {
        Tfar = T2;
      }
    }
    else
    {
      // parallel to an axis - check for intersection at all
      if ((StartIndex[i] < ImStart[i]) || (StartIndex[i] > ImStart[i] + static_cast<int>(ImSize[i]) - 1))
      {
        // no intersection
        start = end = 0;
        return (0);
      }
    }
  }
  sPos = static_cast<int>(Tnear * itk::Math::abs(line[perpdir]) + 0.5);
  ePos = static_cast<int>(Tfar * itk::Math::abs(line[perpdir]) + 0.5);

  // std::cout << Tnear << ' ' << Tfar << std::endl;
  if (Tfar < Tnear) // seems to need some margin
  {
    // in theory, no intersection, but search between them
    bool         intersection = false;
    unsigned int inside = 0; // initialize to avoid warning
    if (Tnear - Tfar < 10)
    {
      //      std::cout << "Searching " << Tnear << ' ' << Tfar << std::endl;
      itkAssertInDebugAndIgnoreInReleaseMacro(ePos >= 0);
      itkAssertInDebugAndIgnoreInReleaseMacro(sPos < static_cast<int>(LineOffsets.size()));
      for (int i = ePos; i <= sPos; ++i)
      {
        if (AllImage.IsInside(StartIndex + LineOffsets[i]))
        {
          inside = i;
          intersection = true;
          break;
        }
      }
    }
    if (intersection)
    {
      //      std::cout << "Found intersection after all :: " << inside << std::endl;
      sPos = ePos = inside;
      itkAssertInDebugAndIgnoreInReleaseMacro(ePos + 1 >= 0);
      itkAssertInDebugAndIgnoreInReleaseMacro(ePos + 1 < static_cast<int>(LineOffsets.size()));
      while (AllImage.IsInside(StartIndex + LineOffsets[ePos + 1]))
      {
        ++ePos;
        itkAssertInDebugAndIgnoreInReleaseMacro(ePos + 1 >= 0);
        itkAssertInDebugAndIgnoreInReleaseMacro(ePos + 1 < static_cast<int>(LineOffsets.size()));
      }
      itkAssertInDebugAndIgnoreInReleaseMacro(sPos - 1 >= 0);
      itkAssertInDebugAndIgnoreInReleaseMacro(sPos - 1 < static_cast<int>(LineOffsets.size()));
      while (AllImage.IsInside(StartIndex + LineOffsets[sPos - 1]))
      {
        --sPos;
        itkAssertInDebugAndIgnoreInReleaseMacro(sPos - 1 >= 0);
        itkAssertInDebugAndIgnoreInReleaseMacro(sPos - 1 < static_cast<int>(LineOffsets.size()));
      }
      start = sPos;
      end = ePos;
    }
    else
    {
      //      std::cout << StartIndex << "No intersection" << std::endl;
      start = end = 0;
      return (0);
    }
  }
  else
  {
    itkAssertInDebugAndIgnoreInReleaseMacro(sPos >= 0);
    itkAssertInDebugAndIgnoreInReleaseMacro(sPos < static_cast<int>(LineOffsets.size()));
    if (AllImage.IsInside(StartIndex + LineOffsets[sPos]))
    {
      for (; sPos > 0;)
      {
        itkAssertInDebugAndIgnoreInReleaseMacro(sPos - 1 >= 0);
        itkAssertInDebugAndIgnoreInReleaseMacro(sPos - 1 < static_cast<int>(LineOffsets.size()));
        if (!AllImage.IsInside(StartIndex + LineOffsets[sPos - 1]))
        {
          break;
        }
        else
        {
          --sPos;
        }
      }
    }
    else
    {
      for (; sPos < static_cast<int>(LineOffsets.size());)
      {
        itkAssertInDebugAndIgnoreInReleaseMacro(sPos >= 0);
        itkAssertInDebugAndIgnoreInReleaseMacro(sPos < static_cast<int>(LineOffsets.size()));
        ++sPos;
        if (!AllImage.IsInside(StartIndex + LineOffsets[sPos]))
        {
          ++sPos;
        }
        else
        {
          break;
        }
      }
    }
    if (AllImage.IsInside(StartIndex + LineOffsets[ePos]))
    {
      for (; ePos < static_cast<int>(LineOffsets.size());)
      {
        itkAssertInDebugAndIgnoreInReleaseMacro(ePos + 1 >= 0);
        itkAssertInDebugAndIgnoreInReleaseMacro(ePos + 1 < static_cast<int>(LineOffsets.size()));
        if (!AllImage.IsInside(StartIndex + LineOffsets[ePos + 1]))
        {
          break;
        }
        else
        {
          ++ePos;
        }
      }
    }
    else
    {
      for (; ePos > 0;)
      {
        --ePos;
        itkAssertInDebugAndIgnoreInReleaseMacro(ePos >= 0);
        itkAssertInDebugAndIgnoreInReleaseMacro(ePos < static_cast<int>(LineOffsets.size()));
        if (!AllImage.IsInside(StartIndex + LineOffsets[ePos]))
        {
          --ePos;
        }
        else
        {
          break;
        }
      }
    }
  }
  start = sPos;
  end = ePos;
  return (1);
}

template <typename TImage, typename TBres>
void
CopyLineToImage(const typename TImage::Pointer            output,
                const typename TImage::IndexType          StartIndex,
                const typename TBres::OffsetArray         LineOffsets,
                std::vector<typename TImage::PixelType> & outbuffer,
                const unsigned int                        start,
                const unsigned int                        end)
{
  unsigned int size = end - start + 1;

  for (unsigned int i = 0; i < size; ++i)
  {
    // itkAssertInDebugAndIgnoreInReleaseMacro(start + i >= 0);
    itkAssertInDebugAndIgnoreInReleaseMacro(start + i < LineOffsets.size());
    output->SetPixel(StartIndex + LineOffsets[start + i], outbuffer[i + 1]); // compat
  }
}

template <typename TInputImage, typename TLine>
typename TInputImage::RegionType
MakeEnlargedFace(const typename TInputImage::ConstPointer itkNotUsed(input),
                 const typename TInputImage::RegionType   AllImage,
                 const TLine                              line)
{
  // the face list calculator strategy fails in multithreaded mode
  // with 1D kernels
  // because it doesn't return faces of the sub-blocks if they don't
  // fall along the edge of the image
  using RegionType = typename TInputImage::RegionType;
  using SizeType = typename TInputImage::SizeType;
  using IndexType = typename TInputImage::IndexType;
  using FaceListType = std::list<RegionType>;
  FaceListType faceList;

  for (unsigned int i = 0; i < TInputImage::ImageDimension; ++i)
  {
    RegionType R1, R2;
    SizeType   S1 = AllImage.GetSize();
    IndexType  I2 = AllImage.GetIndex();

    S1[i] = 1;
    R1 = AllImage;
    R2 = AllImage;

    // the first face will have the same starting index and one
    // dimension removed
    R1.SetSize(S1);

    I2[i] = I2[i] + AllImage.GetSize()[i] - 1;
    R2.SetSize(S1);
    R2.SetIndex(I2);
    faceList.push_back(R1);
    faceList.push_back(R2);
    //    std::cout << R1 << R2 << std::endl;
  }
  typename FaceListType::iterator fit;
  fit = faceList.begin();

  typename TInputImage::RegionType RelevantRegion;
  bool                             foundFace = false;
  float                            MaxComp = NumericTraits<float>::NonpositiveMin();
  unsigned int                     DomDir = 0;
  // std::cout << "------------" << std::endl;
  // figure out the dominant direction of the line
  for (unsigned int i = 0; i < TInputImage::RegionType::ImageDimension; ++i)
  {
    const auto abs_line_elmt_tmp = itk::Math::abs(line[i]);
    if (abs_line_elmt_tmp > MaxComp)
    {
      MaxComp = abs_line_elmt_tmp;
      DomDir = i;
    }
  }

  for (; fit != faceList.end(); ++fit)
  {
    // check whether this face is suitable for parallel sweeping - i.e
    // whether the line is within 45 degrees of the perpendicular
    // Figure out the perpendicular using the region size
    unsigned int FaceDir = 0;
    //    std::cout << "Face " << *fit << std::endl;
    for (unsigned int i = 0; i < TInputImage::RegionType::ImageDimension; ++i)
    {
      if (fit->GetSize()[i] == 1)
      {
        FaceDir = i;
      }
    }
    if (FaceDir == DomDir) // within 1 degree
    {
      // now check whether the line goes inside the image from this face
      if (NeedToDoFace<typename TInputImage::RegionType, TLine>(AllImage, *fit, line))
      {
        //        std::cout << "Using face: " << *fit << line << std::endl;
        RelevantRegion = *fit;
        foundFace = true;
        break;
      }
    }
  }
  if (foundFace)
  {
    // enlarge the region so that sweeping the line across it will
    // cause all pixels to be visited.
    // find the dimension not within the face
    unsigned int NonFaceDim = 0;

    for (unsigned int i = 0; i < TInputImage::RegionType::ImageDimension; ++i)
    {
      if (RelevantRegion.GetSize()[i] == 1)
      {
        NonFaceDim = i;
        break;
      }
    }

    // figure out how much extra each other dimension needs to be extended
    typename TInputImage::SizeType  NewSize = RelevantRegion.GetSize();
    typename TInputImage::IndexType NewStart = RelevantRegion.GetIndex();
    unsigned int                    NonFaceLen = AllImage.GetSize()[NonFaceDim];
    for (unsigned int i = 0; i < TInputImage::RegionType::ImageDimension; ++i)
    {
      if (i != NonFaceDim)
      {
        auto Pad = Math::Ceil<int>(static_cast<float>(NonFaceLen) * line[i] / itk::Math::abs(line[NonFaceDim]));
        if (Pad < 0)
        {
          // just increase the size - no need to change the start
          NewSize[i] += itk::Math::abs(Pad) + 1;
        }
        else
        {
          // change the size and index
          NewSize[i] += Pad + 1;
          NewStart[i] -= Pad + 1;
        }
      }
    }
    RelevantRegion.SetSize(NewSize);
    RelevantRegion.SetIndex(NewStart);
  }
  else
  {
    std::cout << "Line " << line << " doesn't correspond to a face" << std::endl;
  }

  return RelevantRegion;
}

template <typename TImage, typename TBres, typename TLine>
int
FillLineBuffer(typename TImage::ConstPointer             input,
               const typename TImage::IndexType          StartIndex,
               const TLine                               line, // unit vector
               const float                               tol,
               const typename TBres::OffsetArray         LineOffsets,
               const typename TImage::RegionType         AllImage,
               std::vector<typename TImage::PixelType> & inbuffer,
               unsigned int &                            start,
               unsigned int &                            end)
{
  int status = ComputeStartEnd<TImage, TBres, TLine>(StartIndex, line, tol, LineOffsets, AllImage, start, end);
  if (!status)
  {
    return (status);
  }
  unsigned int size = end - start + 1;
  // compat
  for (unsigned int i = 0; i < size; ++i)
  {
    itkAssertInDebugAndIgnoreInReleaseMacro(start + i < LineOffsets.size());
    inbuffer[i + 1] = input->GetPixel(StartIndex + LineOffsets[start + i]);
  }
  return (1);
}

template <typename TLine>
unsigned int
GetLinePixels(const TLine line)
{
  float N = line.GetNorm();
  float correction = 0.0;

  for (unsigned int i = 0; i < TLine::Dimension; ++i)
  {
    const float tt = itk::Math::abs(line[i] / N);
    if (tt > correction)
    {
      correction = tt;
    }
  }

  N *= correction;
  return static_cast<int>(N + 0.5);
}
} // namespace itk

#endif