/*=========================================================================
 *
 *  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 itkVoronoiPartitioningImageFilter_hxx
#define itkVoronoiPartitioningImageFilter_hxx
#include "itkVoronoiPartitioningImageFilter.h"

#include "itkImageRegionIteratorWithIndex.h"

namespace itk
{
template <typename TInputImage, typename TOutputImage>
void
VoronoiPartitioningImageFilter<TInputImage, TOutputImage>::ClassifyDiagram()
{
  PointIdIterator currPit;
  PointIdIterator currPitEnd;
  PointType       currP;
  PointTypeDeque  VertList;
  IndexList       PixelPool;

  this->m_NumberOfBoundary = 0;
  for (int i = 0; i < this->GetNumberOfSeeds(); i++)
  {
    CellAutoPointer currCell;
    this->m_WorkingVD->GetCellId(i, currCell);
    currPitEnd = currCell->PointIdsEnd();
    VertList.clear();
    for (currPit = currCell->PointIdsBegin(); currPit != currPitEnd; ++currPit)
    {
      this->m_WorkingVD->GetPoint((*currPit), &(currP));
      VertList.push_back(currP);
    }

    PixelPool.clear();
    this->GetPixelIndexFromPolygon(VertList, &PixelPool);
    this->m_NumberOfPixels[i] = static_cast<SizeValueType>(PixelPool.size());
    this->m_Label[i] = this->TestHomogeneity(PixelPool);

    // when partitioning the NumberOfBoundary is the number of regions that
    // are not homogeneous
    if (!this->m_Label[i])
    {
      this->m_NumberOfBoundary++;
    }
  }
}

template <typename TInputImage, typename TOutputImage>
void
VoronoiPartitioningImageFilter<TInputImage, TOutputImage>::GenerateAddingSeeds()
{
  EdgeIterator eit;
  auto         eitend = this->m_WorkingVD->EdgeEnd();
  PointType    adds;

  Point<int, 2> seeds;

  // Walk the edges of the diagram, if a seed to either side is
  // no homogeneous, then place a new seed along the middle of the edge
  for (eit = this->m_WorkingVD->EdgeBegin(); eit != eitend; ++eit)
  {
    seeds = this->m_WorkingVD->GetSeedsIDAroundEdge(&*eit);
    // if either seed is not homogeneous
    if ((!this->m_Label[seeds[0]] || !this->m_Label[seeds[1]]) &&
        (this->m_NumberOfPixels[seeds[0]] > this->GetMinRegion()) &&
        (this->m_NumberOfPixels[seeds[1]] > this->GetMinRegion()))
    {
      adds[0] = (eit->m_Left[0] + eit->m_Right[0]) * 0.5;
      adds[1] = (eit->m_Left[1] + eit->m_Right[1]) * 0.5;
      this->m_SeedsToAdded.push_back(adds);
    }
  }
}

template <typename TInputImage, typename TOutputImage>
void
VoronoiPartitioningImageFilter<TInputImage, TOutputImage>::MakeSegmentBoundary()
{
  RegionType region = this->GetInput()->GetRequestedRegion();

  itk::ImageRegionIteratorWithIndex<OutputImageType> oit(this->GetOutput(), region);
  while (!oit.IsAtEnd())
  {
    oit.Set(0);
    ++oit;
  }

  NeighborIdIterator nit;
  NeighborIdIterator nitend;
  for (int i = 0; i < this->GetNumberOfSeeds(); i++)
  {
    nitend = this->m_WorkingVD->NeighborIdsEnd(i);
    for (nit = this->m_WorkingVD->NeighborIdsBegin(i); nit != nitend; ++nit)
    {
      if ((*nit) > i)
      {
        this->drawLine(this->m_WorkingVD->GetSeed(i), this->m_WorkingVD->GetSeed(*nit));
      }
    }
  }
}

template <typename TInputImage, typename TOutputImage>
void
VoronoiPartitioningImageFilter<TInputImage, TOutputImage>::MakeSegmentObject()
{
  RegionType region = this->GetInput()->GetRequestedRegion();

  itk::ImageRegionIteratorWithIndex<OutputImageType> oit(this->GetOutput(), region);
  while (!oit.IsAtEnd())
  {
    oit.Set(0);
    ++oit;
  }
  PointIdIterator currPit;
  PointIdIterator currPitEnd;
  PointType       currP;
  PointTypeDeque  VertList;
  for (int i = 0; i < this->GetNumberOfSeeds(); i++)
  {
    CellAutoPointer currCell;
    this->m_WorkingVD->GetCellId(i, currCell);
    currPitEnd = currCell->PointIdsEnd();
    VertList.clear();
    for (currPit = currCell->PointIdsBegin(); currPit != currPitEnd; ++currPit)
    {
      this->m_WorkingVD->GetPoint((*currPit), &(currP));
      VertList.push_back(currP);
    }
    // Need to fill with an segment identifier
    this->FillPolygon(VertList, static_cast<OutputPixelType>(i));
  }
}

template <typename TInputImage, typename TOutputImage>
bool
VoronoiPartitioningImageFilter<TInputImage, TOutputImage>::TestHomogeneity(IndexList & Plist)
{
  auto          num = static_cast<SizeValueType>(Plist.size());
  SizeValueType i;
  double        getp;
  double        addp = 0;
  double        addpp = 0;

  const InputImageType * inputImage = this->GetInput();

  for (i = 0; i < num; i++)
  {
    getp = (double)(inputImage->GetPixel(Plist[i]));
    addp = addp + getp;
    addpp = addpp + getp * getp;
  }

  double savevar;
  if (num > 1)
  {
    savevar = std::sqrt((addpp - (addp * addp) / static_cast<double>(num)) / (static_cast<double>(num) - 1.0));
  }
  else
  {
    savevar = -1;
  }

  return (savevar >= 0 && std::sqrt(savevar) < m_SigmaThreshold);
}

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

  os << indent << "SigmaThreshold: " << m_SigmaThreshold << std::endl;
}
} // namespace itk

#endif
