// -*- c++ -*-
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkGradientFilter.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.

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

/*----------------------------------------------------------------------------
 Copyright (c) Sandia Corporation
 See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
----------------------------------------------------------------------------*/

#include "vtkGradientFilter.h"

#include "vtkCell.h"
#include "vtkCellData.h"
#include "vtkCellDataToPointData.h"
#include "vtkDataArray.h"
#include "vtkDataSet.h"
#include "vtkIdList.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkUnstructuredGrid.h"

//-----------------------------------------------------------------------------

vtkCxxRevisionMacro(vtkGradientFilter, "1.11");
vtkStandardNewMacro(vtkGradientFilter);

template<class data_type>
void vtkGradientFilterDoComputePointGradients(vtkDataSet *structure,
                                              data_type *scalars,
                                              data_type *gradients);
template<class data_type>
int vtkGradientFilterAddCellContribution(vtkIdType pointId,
                                         double *pointCoord, vtkCell *cell,
                                         data_type *scalars, data_type *g);

template<class data_type>
void vtkGradientFilterDoComputeCellGradients(vtkDataSet *structure,
                                             data_type *scalars,
                                             data_type *gradients);

//-----------------------------------------------------------------------------

vtkGradientFilter::vtkGradientFilter()
{
  this->ResultArrayName = NULL;
  this->FasterApproximation = 0;
  this->SetInputScalars(vtkDataObject::FIELD_ASSOCIATION_POINTS_THEN_CELLS,
                        vtkDataSetAttributes::SCALARS);
}

vtkGradientFilter::~vtkGradientFilter()
{
  this->SetResultArrayName(NULL);
}

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

  os << indent << "Result Array Name:"
     << (this->ResultArrayName ? this->ResultArrayName : "Gradients") << endl;
  os << indent << "Faster Approximation:" << this->FasterApproximation << endl;
}

//-----------------------------------------------------------------------------

void vtkGradientFilter::SetInputScalars(int fieldAssociation, const char *name)
{
  if (   (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_POINTS)
      && (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_CELLS)
      && (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_POINTS_THEN_CELLS) )
    {
    vtkErrorMacro("Input scalars must be associated with points or cells.");
    return;
    }

  this->SetInputArrayToProcess(0, 0, 0, fieldAssociation, name);
}

void vtkGradientFilter::SetInputScalars(int fieldAssociation,
                                        int fieldAttributeType)
{
  if (   (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_POINTS)
      && (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_CELLS)
      && (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_POINTS_THEN_CELLS) )
    {
    vtkErrorMacro("Input scalars must be associated with points or cells.");
    return;
    }

  this->SetInputArrayToProcess(0, 0, 0, fieldAssociation, fieldAttributeType);
}

//-----------------------------------------------------------------------------

int vtkGradientFilter::RequestUpdateExtent(vtkInformation *vtkNotUsed(request),
                                           vtkInformationVector **inputVector,
                                           vtkInformationVector *outputVector)
{
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  // Technically, this code is only correct for pieces extent types.  However,
  // since this class is pretty inefficient for data types that use 3D extents,
  // we'll punt on the ghost levels for them, too.
  int piece, numPieces, ghostLevels;
  
  piece = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
  numPieces = outInfo->Get(
                   vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
  ghostLevels = outInfo->Get(
             vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());
  
  if (numPieces > 1)
    {
    ++ghostLevels;
    }

  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(), piece);
  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(),
              numPieces);
  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
              ghostLevels);
  inInfo->Set(vtkStreamingDemandDrivenPipeline::EXACT_EXTENT(), 1);

  return 1;
}


//-----------------------------------------------------------------------------

static int vtkGradientFilterHasArray(vtkFieldData *fieldData,
                                     vtkDataArray *array)
{
  int numarrays = fieldData->GetNumberOfArrays();
  for (int i = 0; i < numarrays; i++)
    {
    if (array == fieldData->GetArray(i))
      {
      return 1;
      }
    }
  return 0;
}

//-----------------------------------------------------------------------------

int vtkGradientFilter::RequestData(vtkInformation *vtkNotUsed(request),
                                   vtkInformationVector **inputVector,
                                   vtkInformationVector *outputVector)
{
  vtkDebugMacro("RequestData");

  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  vtkDataSet *input
    = vtkDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkDataSet *output
    = vtkDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkDataArray *scalars = this->GetInputArrayToProcess(0, inputVector);

  if (scalars == NULL)
    {
    vtkErrorMacro("No input scalars.");
    return 0;
    }
  if (scalars->GetNumberOfComponents() != 1)
    {
    vtkErrorMacro("Input scalars must have one component.");
    return 0;
    }

  int fieldAssociation;
  if (vtkGradientFilterHasArray(input->GetPointData(), scalars))
    {
    fieldAssociation = vtkDataObject::FIELD_ASSOCIATION_POINTS;
    }
  else if (vtkGradientFilterHasArray(input->GetCellData(), scalars))
    {
    fieldAssociation = vtkDataObject::FIELD_ASSOCIATION_CELLS;
    }
  else
    {
    vtkErrorMacro("Input scalars do not seem to be either point or cell scalars.");
    return 0;
    }

  output->CopyStructure(input);
  output->GetPointData()->PassData(input->GetPointData());
  output->GetCellData()->PassData(input->GetCellData());

  vtkDataArray *gradients
    = vtkDataArray::CreateDataArray(scalars->GetDataType());
  gradients->SetNumberOfComponents(3);
  gradients->SetNumberOfTuples(scalars->GetNumberOfTuples());
  if (this->ResultArrayName)
    {
    gradients->SetName(this->ResultArrayName);
    }
  else
    {
    gradients->SetName("Gradients");
    }

  if (fieldAssociation == vtkDataObject::FIELD_ASSOCIATION_POINTS)
    {
    if (!this->FasterApproximation)
      {
      switch (scalars->GetDataType())
        {
        vtkTemplateMacro(vtkGradientFilterDoComputePointGradients(
                          input,
                          static_cast<VTK_TT *>(scalars->GetVoidPointer(0)),
                          static_cast<VTK_TT *>(gradients->GetVoidPointer(0))));
        }

      output->GetPointData()->AddArray(gradients);
      }
    else // this->FasterApproximation
      {
      // The cell computation is faster and works off of point data anyway.  The
      // faster approximation is to use the cell algorithm and then convert the
      // result to point data.
      vtkDataArray *cellGradients
        = vtkDataArray::CreateDataArray(gradients->GetDataType());
      cellGradients->SetName(gradients->GetName());
      cellGradients->SetNumberOfComponents(3);
      cellGradients->SetNumberOfTuples(input->GetNumberOfCells());

      switch (scalars->GetDataType())
        {
        vtkTemplateMacro(vtkGradientFilterDoComputeCellGradients(
                      input,
                      static_cast<VTK_TT *>(scalars->GetVoidPointer(0)),
                      static_cast<VTK_TT *>(cellGradients->GetVoidPointer(0))));
        }

      // We need to convert cell scalars to points scalars.
      vtkDataSet *dummy = input->NewInstance();
      dummy->CopyStructure(input);
      dummy->GetCellData()->AddArray(cellGradients);

      vtkCellDataToPointData *cd2pd = vtkCellDataToPointData::New();
      cd2pd->SetInput(dummy);
      cd2pd->PassCellDataOff();
      cd2pd->Update();

      // Set the gradients array in the output and cleanup.
      vtkDataArray *pointGradients
        = cd2pd->GetOutput()->GetPointData()->GetArray(gradients->GetName());
      output->GetPointData()->AddArray(pointGradients);
      cd2pd->Delete();
      dummy->Delete();
      cellGradients->Delete();
      }
    }
  else  // fieldAssocation == vtkDataObject::FIELD_ASSOCIATION_CELLS
    {
    // We need to convert cell scalars to points scalars.
    vtkDataSet *dummy = input->NewInstance();
    dummy->CopyStructure(input);
    dummy->GetCellData()->SetScalars(scalars);

    vtkCellDataToPointData *cd2pd = vtkCellDataToPointData::New();
    cd2pd->SetInput(dummy);
    cd2pd->PassCellDataOff();
    cd2pd->Update();
    vtkDataArray *pointScalars
      = cd2pd->GetOutput()->GetPointData()->GetScalars();
    pointScalars->Register(this);
    cd2pd->Delete();
    dummy->Delete();

    switch (pointScalars->GetDataType())
      {
      vtkTemplateMacro(vtkGradientFilterDoComputeCellGradients(
                          input,
                          static_cast<VTK_TT *>(pointScalars->GetVoidPointer(0)),
                          static_cast<VTK_TT *>(gradients->GetVoidPointer(0))));
      }

    output->GetCellData()->AddArray(gradients);
    pointScalars->UnRegister(this);
    }

  gradients->Delete();

  // If necessary, remove a layer of ghost cells.
  int numPieces
    = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
  if (numPieces > 1)
    {
    int ghostLevel = outInfo->Get(
             vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());
    vtkPolyData *pd = vtkPolyData::SafeDownCast(output);
    vtkUnstructuredGrid *ug = vtkUnstructuredGrid::SafeDownCast(output);
    if (pd) pd->RemoveGhostCells(ghostLevel+1);
    if (ug) ug->RemoveGhostCells(ghostLevel+1);
    }

  return 1;
}

//-----------------------------------------------------------------------------

template<class data_type>
void vtkGradientFilterDoComputePointGradients(vtkDataSet *structure,
                                              data_type *scalars,
                                              data_type *gradients)
{
  vtkIdType numpts;
  data_type *g;
  vtkIdType point;
  vtkIdList *currentPoint;
  vtkIdList *cellsOnPoint;

  currentPoint = vtkIdList::New();
  currentPoint->SetNumberOfIds(1);
  cellsOnPoint = vtkIdList::New();

  numpts = structure->GetNumberOfPoints();
  g = gradients;
  for (point = 0; point < numpts; point++)
    {
    g[0] = g[1] = g[2] = 0;
    currentPoint->SetId(0, point);

    double pointcoords[3];
    structure->GetPoint(point, pointcoords);

    // Get all cells touching this point.
    structure->GetCellNeighbors(-1, currentPoint, cellsOnPoint);
    vtkIdType numCellNeighbors = cellsOnPoint->GetNumberOfIds();
    vtkIdType numValidCellNeighbors = 0;

    // Iterate on all cells and find all points connected to current point
    // by an edge.
    for (vtkIdType neighbor = 0; neighbor < numCellNeighbors; neighbor++)
      {
      vtkCell *cell = structure->GetCell(cellsOnPoint->GetId(neighbor));

      numValidCellNeighbors += ::vtkGradientFilterAddCellContribution(
                                          point, pointcoords, cell, scalars, g);
      }

    if (numCellNeighbors > 0)
      {
      g[0] /= numCellNeighbors;
      g[1] /= numCellNeighbors;
      g[2] /= numCellNeighbors;
      }
    g += 3;
    }

  currentPoint->Delete();
  cellsOnPoint->Delete();
}

//-----------------------------------------------------------------------------

template<class data_type>
int vtkGradientFilterAddCellContribution(vtkIdType pointId,
                                         double *pointCoord, vtkCell *cell,
                                         data_type *scalars, data_type *g)
{
  double parametricCoord[3];
  int subId;
  double dummy;
  int numpoints = cell->GetNumberOfPoints();
  double *values = new double[numpoints];
  double derivative[3];
  int i;

  // Watch out for degenerate cells.  They make the derivative calculation
  // fail.
  vtkIdList *pointIds = cell->GetPointIds();
  int timesPointRegistered = 0;
  for (i = 0; i < pointIds->GetNumberOfIds(); i++)
    {
    if (pointId == pointIds->GetId(i)) timesPointRegistered++;
    }
  if (timesPointRegistered != 1)
    {
    // The cell should have the point exactly once.  Not good.
    return 0;
    }

  // Get parametric position of point.
  cell->EvaluatePosition(pointCoord, NULL, subId, parametricCoord,
                         dummy, values/*Really another dummy.*/);

  // Get values of scalars at cell points.
  for (i = 0; i < numpoints; i++)
    {
    values[i] = static_cast<double>(scalars[cell->GetPointId(i)]);
    }

  // Get derivitive of cell at point.
  cell->Derivatives(subId, parametricCoord, values, 1, derivative);

  g[0] += static_cast<data_type>(derivative[0]);
  g[1] += static_cast<data_type>(derivative[1]);
  g[2] += static_cast<data_type>(derivative[2]);

  delete [] values;

  return 1;
}

//-----------------------------------------------------------------------------

template<class data_type>
void vtkGradientFilterDoComputeCellGradients(vtkDataSet *structure,
                                             data_type *scalars,
                                             data_type *gradients)
{
  vtkIdType numcells;
  data_type *g;
  vtkIdType cellid;

  numcells = structure->GetNumberOfCells();
  g = gradients;
  for (cellid = 0; cellid < numcells; cellid++)
    {
    vtkCell *cell = structure->GetCell(cellid);

    int subId;
    double cellCenter[3];
    subId = cell->GetParametricCenter(cellCenter);

    int numpoints = cell->GetNumberOfPoints();
    double derivative[3];
    double *values = new double[numpoints];
    for (int i = 0; i < numpoints; i++)
      {
      values[i] = static_cast<double>(scalars[cell->GetPointId(i)]);
      }

    cell->Derivatives(subId, cellCenter, values, 1, derivative);
    delete[] values;
    g[0] = static_cast<data_type>(derivative[0]);
    g[1] = static_cast<data_type>(derivative[1]);
    g[2] = static_cast<data_type>(derivative[2]);
    g += 3;
    }
}
