/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkCellDataToPointData.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

  This software is distributed WITHOUT ANY WARRANTY; without even
  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  PURPOSE.  See the above copyright notice for more information.

  =========================================================================*/
#include "vtkCellDataToPointData.h"

#include "vtkCellData.h"
#include "vtkCell.h"
#include "vtkDataSet.h"
#include "vtkIdList.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkSmartPointer.h"
#include "vtkStructuredGrid.h"
#include "vtkUniformGrid.h"
#include "vtkUnsignedIntArray.h"

#include <algorithm>
#include <functional>
#include <set>

#define VTK_MAX_CELLS_PER_POINT 4096

vtkStandardNewMacro(vtkCellDataToPointData);

namespace
{
//----------------------------------------------------------------------------
// Helper template function that implement the major part of the algorighm
// which will be expanded by the vtkTemplateMacro. The template function is
// provided so that coverage test can cover this function.
  template <typename T>
  void __spread (vtkDataSet* const src, vtkUnsignedIntArray* const num,
                 vtkDataArray* const srcarray, vtkDataArray* const dstarray,
                 vtkIdType ncells, vtkIdType npoints, vtkIdType ncomps,
                 int highestCellDimension, int contributingCellOption)
  {
    T const* const srcptr = static_cast<T const*>(srcarray->GetVoidPointer(0));
    T      * const dstptr = static_cast<T      *>(dstarray->GetVoidPointer(0));

    // zero initialization
    std::fill_n(dstptr, npoints*ncomps, T(0));

    // accumulate
    if (contributingCellOption != vtkCellDataToPointData::Patch)
    {
      T const* srcbeg = srcptr;
      for (vtkIdType cid = 0; cid < ncells; ++cid, srcbeg += ncomps)
      {
        vtkCell* cell = src->GetCell(cid);
        if (cell->GetCellDimension() >= highestCellDimension)
        {
          vtkIdList* pids = cell->GetPointIds();
          for (vtkIdType i = 0, I = pids->GetNumberOfIds(); i < I; ++i)
          {
            T* const dstbeg = dstptr + pids->GetId(i)*ncomps;
            // accumulate cell data to point data <==> point_data += cell_data
            std::transform(srcbeg,srcbeg+ncomps,dstbeg,dstbeg,std::plus<T>());
          }
        }
      }
      // average
      T* dstbeg = dstptr;
      for (vtkIdType pid = 0; pid < npoints; ++pid, dstbeg += ncomps)
      {
        // guard against divide by zero
        if (unsigned int const denom = num->GetValue(pid))
        {
          // divide point data by the number of cells using it <==>
          // point_data /= denum
          std::transform(dstbeg, dstbeg+ncomps, dstbeg,
                         std::bind(std::divides<T>(), std::placeholders::_1, denom));
        }
      }
    }
    else
    { // compute over cell patches
      vtkNew<vtkIdList> cellsOnPoint;
      std::vector<T> data(4*ncomps);
      for (vtkIdType pid = 0; pid < npoints; ++pid)
      {
        std::fill(data.begin(), data.end(), 0);
        T numPointCells[4] = {0, 0, 0, 0};
        // Get all cells touching this point.
        src->GetPointCells(pid, cellsOnPoint);
        vtkIdType numPatchCells = cellsOnPoint->GetNumberOfIds();
        for (vtkIdType pc = 0; pc  < numPatchCells; pc++)
        {
          vtkIdType cellId = cellsOnPoint->GetId(pc);
          int cellDimension = src->GetCell(cellId)->GetCellDimension();
          numPointCells[cellDimension] += 1;
          for (int comp=0;comp<ncomps;comp++)
          {
            data[comp+ncomps*cellDimension] += srcptr[comp+cellId*ncomps];
          }
        }
        for (int dimension=3;dimension>=0;dimension--)
        {
          if (numPointCells[dimension])
          {
            for (int comp=0;comp<ncomps;comp++)
            {
              dstptr[comp+pid*ncomps] = data[comp+dimension*ncomps] / numPointCells[dimension];
            }
            break;
          }
        }
      }
    }
  }
} // end anonymous namespace

class vtkCellDataToPointData::Internals
{
  public:
    std::set<std::string> CellDataArrays;

    // Special traversal algorithm for vtkUniformGrid and vtkRectilinearGrid to support blanking
    // points will not have more than 8 cells for either of these data sets
    template <typename T>
    int InterpolatePointDataWithMask(vtkCellDataToPointData* filter, T *input, vtkDataSet *output)
    {
      vtkNew<vtkIdList> allCellIds;
      allCellIds->Allocate(8);
      vtkNew<vtkIdList> cellIds;
      cellIds->Allocate(8);

      vtkIdType numPts = input->GetNumberOfPoints();

      vtkCellData *inputInCD = input->GetCellData();
      vtkCellData *inCD;
      vtkPointData *outPD = output->GetPointData();

      if (!filter->GetProcessAllArrays())
      {
        inCD = vtkCellData::New();

        for (const auto &name : this->CellDataArrays)
        {
          vtkAbstractArray *arr = inputInCD->GetAbstractArray(name.c_str());
          if (arr == nullptr)
          {
            vtkWarningWithObjectMacro(filter, "cell data array name not found.");
            continue;
          }
          inCD->AddArray(arr);
        }
      }
      else
      {
        inCD = inputInCD;
      }

      outPD->InterpolateAllocate(inCD,numPts);

      double weights[8];

      int abort = 0;
      vtkIdType progressInterval = numPts / 20 + 1;
      for (vtkIdType ptId = 0; ptId < numPts && !abort; ptId++)
      {
        if ( !(ptId % progressInterval) )
        {
          filter->UpdateProgress(static_cast<double>(ptId)/numPts);
          abort = filter->GetAbortExecute();
        }
        input->GetPointCells(ptId, allCellIds);
        cellIds->Reset();
        // Only consider cells that are not masked:
        for (vtkIdType cId = 0; cId < allCellIds->GetNumberOfIds(); ++cId)
        {
          vtkIdType curCell = allCellIds->GetId(cId);
          if (input->IsCellVisible(curCell))
          {
            cellIds->InsertNextId(curCell);
          }
        }

        vtkIdType numCells = cellIds->GetNumberOfIds();

        if ( numCells > 0 )
        {
          double weight = 1.0 / numCells;
          for (vtkIdType cellId = 0; cellId < numCells; cellId++)
          {
            weights[cellId] = weight;
          }
          outPD->InterpolatePoint(inCD, ptId, cellIds, weights);
        }
        else
        {
          outPD->NullPoint(ptId);
        }
      }

      if(!filter->GetProcessAllArrays())
      {
        inCD->Delete();
      }

      return 1;
    }
};

//----------------------------------------------------------------------------
// Instantiate object so that cell data is not passed to output.
vtkCellDataToPointData::vtkCellDataToPointData()
{
  this->PassCellData = 0;
  this->ContributingCellOption = vtkCellDataToPointData::All;
  this->ProcessAllArrays = true;
  this->Implementation = new Internals();
}

//----------------------------------------------------------------------------
vtkCellDataToPointData::~vtkCellDataToPointData()
{
  delete this->Implementation;
}

//----------------------------------------------------------------------------
void vtkCellDataToPointData::AddCellDataArray(const char *name)
{
  if (!name)
  {
    vtkErrorMacro("name cannot be null.");
    return;
  }

  this->Implementation->CellDataArrays.insert(std::string(name));
  this->Modified();
}

//----------------------------------------------------------------------------
void vtkCellDataToPointData::RemoveCellDataArray(const char *name)
{
  if (!name)
  {
    vtkErrorMacro("name cannot be null.");
    return;
  }

  this->Implementation->CellDataArrays.erase(name);
  this->Modified();
}

//----------------------------------------------------------------------------
void vtkCellDataToPointData::ClearCellDataArrays()
{
  if (!this->Implementation->CellDataArrays.empty())
  {
    this->Modified();
  }
  this->Implementation->CellDataArrays.clear();
}

//----------------------------------------------------------------------------
int vtkCellDataToPointData::RequestData(
  vtkInformation*,
  vtkInformationVector** inputVector,
  vtkInformationVector* outputVector)
{
  vtkInformation* info = outputVector->GetInformationObject(0);
  vtkDataSet *output = vtkDataSet::SafeDownCast(
    info->Get(vtkDataObject::DATA_OBJECT()));

  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkDataSet *input = vtkDataSet::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkDebugMacro(<<"Mapping cell data to point data");

  // Special traversal algorithm for unstructured grid
  if (input->IsA("vtkUnstructuredGrid") || input->IsA("vtkPolyData"))
  {
    return this->RequestDataForUnstructuredData(nullptr, inputVector, outputVector);
  }

  vtkDebugMacro(<<"Mapping cell data to point data");

  // First, copy the input to the output as a starting point
  output->CopyStructure( input );

  // Pass the point data first. The fields and attributes
  // which also exist in the cell data of the input will
  // be over-written during CopyAllocate
  output->GetPointData()->CopyGlobalIdsOff();
  output->GetPointData()->PassData(input->GetPointData());
  output->GetPointData()->CopyFieldOff(vtkDataSetAttributes::GhostArrayName());

  if (input->GetNumberOfPoints() < 1)
  {
    vtkDebugMacro(<<"No input point data!");
    return 1;
  }

  // Do the interpolation, taking care of masked cells if needed.
  vtkStructuredGrid *sGrid = vtkStructuredGrid::SafeDownCast(input);
  vtkUniformGrid *uniformGrid = vtkUniformGrid::SafeDownCast(input);
  int result;
  if (sGrid && sGrid->HasAnyBlankCells())
  {
    result = this->Implementation->InterpolatePointDataWithMask(this, sGrid, output);
  }
  else if (uniformGrid && uniformGrid->HasAnyBlankCells())
  {
    result = this->Implementation->InterpolatePointDataWithMask(this, uniformGrid, output);
  }
  else
  {
    result = this->InterpolatePointData(input, output);
  }

  if (result == 0)
  {
    return 0;
  }

  if (!this->PassCellData)
  {
    output->GetCellData()->CopyAllOff();
    output->GetCellData()->CopyFieldOn(vtkDataSetAttributes::GhostArrayName());
  }
  output->GetCellData()->PassData(input->GetCellData());

  return 1;
}

//----------------------------------------------------------------------------
void vtkCellDataToPointData::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);

  os << indent << "PassCellData: " << (this->PassCellData ? "On\n" : "Off\n");
  os << indent << "ContributingCellOption: " << this->ContributingCellOption << endl;
}

//----------------------------------------------------------------------------
int vtkCellDataToPointData::RequestDataForUnstructuredData
(vtkInformation*,
 vtkInformationVector** inputVector,
 vtkInformationVector* outputVector)
{
  vtkDataSet* const src = vtkDataSet::SafeDownCast(
    inputVector[0]->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT()));
  vtkDataSet* const dst = vtkDataSet::SafeDownCast(
    outputVector->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT()));

  vtkIdType const ncells  = src->GetNumberOfCells ();
  vtkIdType const npoints = src->GetNumberOfPoints();
  if (ncells < 1 || npoints < 1)
  {
    vtkDebugMacro(<<"No input data!");
    return 1;
  }

  // count the number of cells associated with each point. if we are doing patches
  // though we will do that later on.
  vtkSmartPointer<vtkUnsignedIntArray> num;
  int highestCellDimension = 0;
  if (this->ContributingCellOption != vtkCellDataToPointData::Patch)
  {
    num = vtkSmartPointer<vtkUnsignedIntArray>::New();
    num->SetNumberOfComponents(1);
    num->SetNumberOfTuples(npoints);
    std::fill_n(num->GetPointer(0), npoints, 0u);
    if (this->ContributingCellOption == vtkCellDataToPointData::DataSetMax)
    {
      int maxDimension = src->IsA("vtkPolyData") == 1 ? 2 : 3;
      for (vtkIdType i=0;i<src->GetNumberOfCells();i++)
      {
        int dim = src->GetCell(i)->GetCellDimension();
        if (dim > highestCellDimension)
        {
          highestCellDimension = dim;
          if (highestCellDimension == maxDimension)
          {
            break;
          }
        }
      }
    }
    vtkNew<vtkIdList> pids;
    for (vtkIdType cid = 0; cid < ncells; ++cid)
    {
      if (src->GetCell(cid)->GetCellDimension() >= highestCellDimension)
      {
        src->GetCellPoints(cid, pids);
        for (vtkIdType i = 0, I = pids->GetNumberOfIds(); i < I; ++i)
        {
          vtkIdType const pid = pids->GetId(i);
          num->SetValue(pid, num->GetValue(pid)+1);
        }
      }
    }
  }

  // First, copy the input to the output as a starting point
  dst->CopyStructure(src);
  vtkPointData* const opd = dst->GetPointData();

  // Pass the point data first. The fields and attributes
  // which also exist in the cell data of the input will
  // be over-written during CopyAllocate
  opd->CopyGlobalIdsOff();
  opd->PassData(src->GetPointData());
  opd->CopyFieldOff(vtkDataSetAttributes::GhostArrayName());

  // Copy all existing cell fields into a temporary cell data array,
  // unless the SelectCellDataArrays option is active
  vtkSmartPointer<vtkCellData> processedCellData;
  if (!this->ProcessAllArrays)
  {
    processedCellData = vtkSmartPointer<vtkCellData>::New();

    vtkCellData *processedCellDataTemp = src->GetCellData();
    for (const auto &name : this->Implementation->CellDataArrays)
    {
      vtkAbstractArray *arr = processedCellDataTemp->GetAbstractArray(name.c_str());
      if (arr == nullptr)
      {
        vtkWarningMacro("cell data array name not found.");
        continue;
      }
      processedCellData->AddArray(arr);
    }
  }
  else
  {
    processedCellData = vtkSmartPointer<vtkCellData>(src->GetCellData());
  }

  // Remove all fields that are not a data array.
  for (vtkIdType fid = processedCellData->GetNumberOfArrays(); fid--;)
  {
    if (!vtkDataArray::FastDownCast(processedCellData->GetAbstractArray(fid)))
    {
      processedCellData->RemoveArray(fid);
    }
  }

  // Cell field list constructed from the filtered cell data array
  vtkDataSetAttributes::FieldList cfl(1);
  cfl.InitializeFieldList(processedCellData);
  opd->InterpolateAllocate(cfl, npoints, npoints);

  const auto nfields = processedCellData->GetNumberOfArrays();
  int fid = 0;
  auto f = [this, &fid, nfields, npoints, src, num, ncells, highestCellDimension](
             vtkAbstractArray* aa_srcarray, vtkAbstractArray* aa_dstarray) {
    // update progress and check for an abort request.
    this->UpdateProgress((fid + 1.0) / nfields);
    ++fid;

    if (this->GetAbortExecute())
    {
      return;
    }

    vtkDataArray* const srcarray = vtkDataArray::FastDownCast(aa_srcarray);
    vtkDataArray* const dstarray = vtkDataArray::FastDownCast(aa_dstarray);
    if (srcarray && dstarray)
    {
      dstarray->SetNumberOfTuples(npoints);
      vtkIdType const ncomps = srcarray->GetNumberOfComponents();
      switch (srcarray->GetDataType())
      {
        vtkTemplateMacro(__spread<VTK_TT>(src, num, srcarray, dstarray, ncells, npoints, ncomps,
          highestCellDimension, this->ContributingCellOption));
      }
    }
  };

  if (processedCellData != nullptr && dst->GetPointData() != nullptr)
  {
    cfl.TransformData(0, processedCellData, dst->GetPointData(), f);
  }

  if (!this->PassCellData)
  {
    dst->GetCellData()->CopyAllOff();
    dst->GetCellData()->CopyFieldOn(vtkDataSetAttributes::GhostArrayName());
  }
  dst->GetCellData()->PassData(src->GetCellData());

  return 1;
}

int vtkCellDataToPointData::InterpolatePointData(vtkDataSet *input, vtkDataSet *output)
{
  vtkNew<vtkIdList> cellIds;
  cellIds->Allocate(VTK_MAX_CELLS_PER_POINT);

  vtkIdType numPts = input->GetNumberOfPoints();

  vtkCellData *inputInCD = input->GetCellData();
  vtkCellData *inCD;
  vtkPointData *outPD = output->GetPointData();

  if (!this->ProcessAllArrays)
  {
    inCD = vtkCellData::New();

    for (const auto &name : this->Implementation->CellDataArrays)
    {
      vtkAbstractArray *arr = inputInCD->GetAbstractArray(name.c_str());
      if (arr == nullptr)
      {
        vtkWarningMacro("cell data array name not found.");
        continue;
      }
      inCD->AddArray(arr);
    }
  }
  else
  {
    inCD = inputInCD;
  }

  outPD->InterpolateAllocate(inCD,numPts);

  double weights[VTK_MAX_CELLS_PER_POINT];

  int abort = 0;
  vtkIdType progressInterval = numPts / 20 + 1;
  for (vtkIdType ptId = 0; ptId < numPts && !abort; ptId++)
  {
    if ( !(ptId % progressInterval) )
    {
      this->UpdateProgress(static_cast<double>(ptId)/numPts);
      abort = GetAbortExecute();
    }

    input->GetPointCells(ptId, cellIds);
    vtkIdType numCells = cellIds->GetNumberOfIds();

    if ( numCells > 0 && numCells < VTK_MAX_CELLS_PER_POINT )
    {
      double weight = 1.0 / numCells;
      for (vtkIdType cellId = 0; cellId < numCells; cellId++)
      {
        weights[cellId] = weight;
      }
      outPD->InterpolatePoint(inCD, ptId, cellIds, weights);
    }
    else
    {
      outPD->NullPoint(ptId);
    }
  }

  if (!this->ProcessAllArrays)
  {
    inCD->Delete();
  }

  return 1;
}
