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

  Program:   Visualization Toolkit
  Module:    vtkExtractSelectedThresholds.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 "vtkExtractSelectedThresholds.h"

#include "vtkArrayDispatch.h"
#include "vtkCell.h"
#include "vtkCellData.h"
#include "vtkDataArrayRange.h"
#include "vtkDataSet.h"
#include "vtkDoubleArray.h"
#include "vtkIdList.h"
#include "vtkIdTypeArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkSelection.h"
#include "vtkSelectionNode.h"
#include "vtkSignedCharArray.h"
#include "vtkSmartPointer.h"
#include "vtkTable.h"
#include "vtkThreshold.h"
#include "vtkUnstructuredGrid.h"

#include <algorithm>

vtkStandardNewMacro(vtkExtractSelectedThresholds);

//----------------------------------------------------------------------------
vtkExtractSelectedThresholds::vtkExtractSelectedThresholds()
{
  this->SetNumberOfInputPorts(2);
}

//----------------------------------------------------------------------------
vtkExtractSelectedThresholds::~vtkExtractSelectedThresholds() = default;

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

  vtkDataObject* inputDO = vtkDataObject::GetData(inInfo);

  // verify the input, selection and output
  if (!selInfo)
  {
    // When not given a selection, quietly select nothing.
    return 1;
  }

  vtkSelection* sel = vtkSelection::GetData(selInfo);
  vtkSelectionNode* node = nullptr;
  if (sel->GetNumberOfNodes() == 1)
  {
    node = sel->GetNode(0);
  }
  if (!node)
  {
    vtkErrorMacro("Selection must have a single node.");
    return 1;
  }
  if (!node->GetProperties()->Has(vtkSelectionNode::CONTENT_TYPE()) ||
    node->GetProperties()->Get(vtkSelectionNode::CONTENT_TYPE()) != vtkSelectionNode::THRESHOLDS)
  {
    vtkErrorMacro("Missing or invalid CONTENT_TYPE.");
    return 1;
  }

  if (vtkDataSet* input = vtkDataSet::SafeDownCast(inputDO))
  {
    if (input->GetNumberOfCells() == 0 && input->GetNumberOfPoints() == 0)
    {
      // empty input, nothing to do..
      return 1;
    }

    vtkDataSet* output = vtkDataSet::GetData(outInfo);
    vtkDebugMacro(<< "Extracting from dataset");

    int thresholdByPointVals = 0;
    int fieldType = vtkSelectionNode::CELL;
    if (node->GetProperties()->Has(vtkSelectionNode::FIELD_TYPE()))
    {
      fieldType = node->GetProperties()->Get(vtkSelectionNode::FIELD_TYPE());
      if (fieldType == vtkSelectionNode::POINT)
      {
        if (node->GetProperties()->Has(vtkSelectionNode::CONTAINING_CELLS()))
        {
          thresholdByPointVals = node->GetProperties()->Get(vtkSelectionNode::CONTAINING_CELLS());
        }
      }
    }

    if (thresholdByPointVals || fieldType == vtkSelectionNode::CELL)
    {
      return this->ExtractCells(node, input, output, thresholdByPointVals);
    }
    if (fieldType == vtkSelectionNode::POINT)
    {
      return this->ExtractPoints(node, input, output);
    }
  }
  else if (vtkTable* inputTable = vtkTable::SafeDownCast(inputDO))
  {
    if (inputTable->GetNumberOfRows() == 0)
    {
      return 1;
    }
    vtkTable* output = vtkTable::GetData(outInfo);
    return this->ExtractRows(node, inputTable, output);
  }

  return 0;
}

//----------------------------------------------------------------------------
int vtkExtractSelectedThresholds::ExtractCells(
  vtkSelectionNode* sel, vtkDataSet* input, vtkDataSet* output, int usePointScalars)
{
  // find the values to threshold within
  vtkDataArray* lims = vtkArrayDownCast<vtkDataArray>(sel->GetSelectionList());
  if (lims == nullptr)
  {
    vtkErrorMacro(<< "No values to threshold with");
    return 1;
  }

  // find out what array we are supposed to threshold in
  vtkDataArray* inScalars = nullptr;
  bool use_ids = false;
  if (usePointScalars)
  {
    if (sel->GetSelectionList()->GetName())
    {
      if (strcmp(sel->GetSelectionList()->GetName(), "vtkGlobalIds") == 0)
      {
        inScalars = input->GetPointData()->GetGlobalIds();
      }
      else if (strcmp(sel->GetSelectionList()->GetName(), "vtkIndices") == 0)
      {
        use_ids = true;
      }
      else
      {
        inScalars = input->GetPointData()->GetArray(sel->GetSelectionList()->GetName());
      }
    }
    else
    {
      inScalars = input->GetPointData()->GetScalars();
    }
  }
  else
  {
    if (sel->GetSelectionList()->GetName())
    {
      if (strcmp(sel->GetSelectionList()->GetName(), "vtkGlobalIds") == 0)
      {
        inScalars = input->GetCellData()->GetGlobalIds();
      }
      else if (strcmp(sel->GetSelectionList()->GetName(), "vtkIndices") == 0)
      {
        use_ids = true;
      }
      else
      {
        inScalars = input->GetCellData()->GetArray(sel->GetSelectionList()->GetName());
      }
    }
    else
    {
      inScalars = input->GetCellData()->GetScalars();
    }
  }
  if (inScalars == nullptr && !use_ids)
  {
    vtkErrorMacro("Could not figure out what array to threshold in.");
    return 1;
  }

  int inverse = 0;
  if (sel->GetProperties()->Has(vtkSelectionNode::INVERSE()))
  {
    inverse = sel->GetProperties()->Get(vtkSelectionNode::INVERSE());
  }

  int passThrough = 0;
  if (this->PreserveTopology)
  {
    passThrough = 1;
  }

  int comp_no = 0;
  if (sel->GetProperties()->Has(vtkSelectionNode::COMPONENT_NUMBER()))
  {
    comp_no = sel->GetProperties()->Get(vtkSelectionNode::COMPONENT_NUMBER());
  }

  vtkIdType cellId, newCellId;
  vtkIdList *cellPts, *pointMap = nullptr;
  vtkIdList* newCellPts = nullptr;
  vtkCell* cell = nullptr;
  vtkPoints* newPoints = nullptr;
  vtkIdType i, ptId, newId, numPts, numCells;
  vtkIdType numCellPts;
  double x[3];

  vtkPointData *pd = input->GetPointData(), *outPD = output->GetPointData();
  vtkCellData *cd = input->GetCellData(), *outCD = output->GetCellData();
  int keepCell;

  outPD->CopyGlobalIdsOn();
  outPD->CopyAllocate(pd);
  outCD->CopyGlobalIdsOn();
  outCD->CopyAllocate(cd);

  numPts = input->GetNumberOfPoints();
  numCells = input->GetNumberOfCells();

  vtkDataSet* outputDS = output;
  vtkSignedCharArray* pointInArray = nullptr;
  vtkSignedCharArray* cellInArray = nullptr;

  vtkUnstructuredGrid* outputUG = nullptr;
  vtkIdTypeArray* originalCellIds = nullptr;
  vtkIdTypeArray* originalPointIds = nullptr;

  signed char flag = inverse ? 1 : -1;

  if (passThrough)
  {
    outputDS->ShallowCopy(input);

    pointInArray = vtkSignedCharArray::New();
    pointInArray->SetNumberOfComponents(1);
    pointInArray->SetNumberOfTuples(numPts);
    for (i = 0; i < numPts; i++)
    {
      pointInArray->SetValue(i, flag);
    }
    pointInArray->SetName("vtkInsidedness");
    outPD->AddArray(pointInArray);
    outPD->SetScalars(pointInArray);

    cellInArray = vtkSignedCharArray::New();
    cellInArray->SetNumberOfComponents(1);
    cellInArray->SetNumberOfTuples(numCells);
    for (i = 0; i < numCells; i++)
    {
      cellInArray->SetValue(i, flag);
    }
    cellInArray->SetName("vtkInsidedness");
    outCD->AddArray(cellInArray);
    outCD->SetScalars(cellInArray);
  }
  else
  {
    outputUG = vtkUnstructuredGrid::SafeDownCast(output);
    outputUG->Allocate(input->GetNumberOfCells());
    newPoints = vtkPoints::New();
    newPoints->Allocate(numPts);

    pointMap = vtkIdList::New(); // maps old point ids into new
    pointMap->SetNumberOfIds(numPts);
    for (i = 0; i < numPts; i++)
    {
      pointMap->SetId(i, -1);
    }

    newCellPts = vtkIdList::New();

    originalCellIds = vtkIdTypeArray::New();
    originalCellIds->SetName("vtkOriginalCellIds");
    originalCellIds->SetNumberOfComponents(1);
    outCD->AddArray(originalCellIds);

    originalPointIds = vtkIdTypeArray::New();
    originalPointIds->SetName("vtkOriginalPointIds");
    originalPointIds->SetNumberOfComponents(1);
    outPD->AddArray(originalPointIds);
    originalPointIds->Delete();
  }

  flag = -flag;

  // Check that the scalars of each cell satisfy the threshold criterion
  for (cellId = 0; cellId < input->GetNumberOfCells(); cellId++)
  {
    cell = input->GetCell(cellId);
    cellPts = cell->GetPointIds();
    numCellPts = cell->GetNumberOfPoints();

    // BUG: This code misses the case where the threshold is contained
    // completely within the cell but none of its points are inside
    // the range.  Consider as an example the threshold range [1, 2]
    // with a cell [0, 3].
    if (usePointScalars)
    {
      keepCell = 0;
      int totalAbove = 0;
      int totalBelow = 0;
      for (i = 0; (i < numCellPts) && (passThrough || !keepCell); i++)
      {
        int above = 0;
        int below = 0;
        ptId = cellPts->GetId(i);
        int inside = this->EvaluateValue(inScalars, comp_no, ptId, lims, &above, &below, nullptr);
        totalAbove += above;
        totalBelow += below;
        // Have we detected a cell that straddles the threshold?
        if ((!inside) && (totalAbove && totalBelow))
        {
          inside = 1;
        }
        if (passThrough && (inside ^ inverse))
        {
          pointInArray->SetValue(ptId, flag);
          cellInArray->SetValue(cellId, flag);
        }
        keepCell |= inside;
      }
    }
    else // use cell scalars
    {
      keepCell = this->EvaluateValue(inScalars, comp_no, cellId, lims);
      if (passThrough && (keepCell ^ inverse))
      {
        cellInArray->SetValue(cellId, flag);
      }
    }

    if (!passThrough && (numCellPts > 0) && (keepCell + inverse == 1)) // Poor man's XOR
    {
      // satisfied thresholding (also non-empty cell, i.e. not VTK_EMPTY_CELL)
      originalCellIds->InsertNextValue(cellId);

      for (i = 0; i < numCellPts; i++)
      {
        ptId = cellPts->GetId(i);
        if ((newId = pointMap->GetId(ptId)) < 0)
        {
          input->GetPoint(ptId, x);
          newId = newPoints->InsertNextPoint(x);
          pointMap->SetId(ptId, newId);
          outPD->CopyData(pd, ptId, newId);
          originalPointIds->InsertNextValue(ptId);
        }
        newCellPts->InsertId(i, newId);
      }
      newCellId = outputUG->InsertNextCell(cell->GetCellType(), newCellPts);
      outCD->CopyData(cd, cellId, newCellId);
      newCellPts->Reset();
    } // satisfied thresholding
  }   // for all cells

  // now clean up / update ourselves
  if (passThrough)
  {
    pointInArray->Delete();
    cellInArray->Delete();
  }
  else
  {
    outputUG->SetPoints(newPoints);
    newPoints->Delete();
    pointMap->Delete();
    newCellPts->Delete();
    originalCellIds->Delete();
  }

  output->Squeeze();

  return 1;
}

//----------------------------------------------------------------------------
int vtkExtractSelectedThresholds::ExtractPoints(
  vtkSelectionNode* sel, vtkDataSet* input, vtkDataSet* output)
{
  // find the values to threshold within
  vtkDataArray* lims = vtkArrayDownCast<vtkDataArray>(sel->GetSelectionList());
  if (lims == nullptr)
  {
    vtkErrorMacro(<< "No values to threshold with");
    return 1;
  }

  // find out what array we are supposed to threshold in
  vtkDataArray* inScalars = nullptr;
  bool use_ids = false;
  if (sel->GetSelectionList()->GetName())
  {
    if (strcmp(sel->GetSelectionList()->GetName(), "vtkGlobalIds") == 0)
    {
      inScalars = input->GetPointData()->GetGlobalIds();
    }
    else if (strcmp(sel->GetSelectionList()->GetName(), "vtkIndices") == 0)
    {
      use_ids = true;
    }
    else
    {
      inScalars = input->GetPointData()->GetArray(sel->GetSelectionList()->GetName());
    }
  }
  else
  {
    inScalars = input->GetPointData()->GetScalars();
  }
  if (inScalars == nullptr && !use_ids)
  {
    vtkErrorMacro("Could not figure out what array to threshold in.");
    return 1;
  }

  int inverse = 0;
  if (sel->GetProperties()->Has(vtkSelectionNode::INVERSE()))
  {
    inverse = sel->GetProperties()->Get(vtkSelectionNode::INVERSE());
  }

  int passThrough = 0;
  if (this->PreserveTopology)
  {
    passThrough = 1;
  }

  int comp_no = 0;
  if (sel->GetProperties()->Has(vtkSelectionNode::COMPONENT_NUMBER()))
  {
    comp_no = sel->GetProperties()->Get(vtkSelectionNode::COMPONENT_NUMBER());
  }

  vtkIdType numPts = input->GetNumberOfPoints();
  vtkPointData* inputPD = input->GetPointData();
  vtkPointData* outPD = output->GetPointData();

  vtkDataSet* outputDS = output;
  vtkSignedCharArray* pointInArray = nullptr;

  vtkUnstructuredGrid* outputUG = nullptr;
  vtkPoints* newPts = vtkPoints::New();

  vtkIdTypeArray* originalPointIds = nullptr;

  signed char flag = inverse ? 1 : -1;

  if (passThrough)
  {
    outputDS->ShallowCopy(input);

    pointInArray = vtkSignedCharArray::New();
    pointInArray->SetNumberOfComponents(1);
    pointInArray->SetNumberOfTuples(numPts);
    for (vtkIdType i = 0; i < numPts; i++)
    {
      pointInArray->SetValue(i, flag);
    }
    pointInArray->SetName("vtkInsidedness");
    outPD->AddArray(pointInArray);
    outPD->SetScalars(pointInArray);
  }
  else
  {
    outputUG = vtkUnstructuredGrid::SafeDownCast(output);
    outputUG->Allocate(numPts);

    newPts->Allocate(numPts);
    outputUG->SetPoints(newPts);

    outPD->CopyGlobalIdsOn();
    outPD->CopyAllocate(inputPD);

    originalPointIds = vtkIdTypeArray::New();
    originalPointIds->SetNumberOfComponents(1);
    originalPointIds->SetName("vtkOriginalPointIds");
    outPD->AddArray(originalPointIds);
    originalPointIds->Delete();
  }

  flag = -flag;

  vtkIdType outPtCnt = 0;
  for (vtkIdType ptId = 0; ptId < numPts; ptId++)
  {
    int keepPoint = this->EvaluateValue(inScalars, comp_no, ptId, lims);
    if (keepPoint ^ inverse)
    {
      if (passThrough)
      {
        pointInArray->SetValue(ptId, flag);
      }
      else
      {
        double X[4];
        input->GetPoint(ptId, X);
        newPts->InsertNextPoint(X);
        outPD->CopyData(inputPD, ptId, outPtCnt);
        originalPointIds->InsertNextValue(ptId);
        outputUG->InsertNextCell(VTK_VERTEX, 1, &outPtCnt);
        outPtCnt++;
      }
    }
  }

  if (passThrough)
  {
    pointInArray->Delete();
  }
  newPts->Delete();
  output->Squeeze();
  return 1;
}

//----------------------------------------------------------------------------
int vtkExtractSelectedThresholds::ExtractRows(
  vtkSelectionNode* sel, vtkTable* input, vtkTable* output)
{
  // find the values to threshold within
  vtkDataArray* lims = vtkArrayDownCast<vtkDataArray>(sel->GetSelectionList());
  if (lims == nullptr)
  {
    vtkErrorMacro(<< "No values to threshold with");
    return 1;
  }

  // Determine the array to threshold.
  vtkDataArray* inScalars = nullptr;
  bool use_ids = false;
  if (sel->GetSelectionList()->GetName())
  {
    if (strcmp(sel->GetSelectionList()->GetName(), "vtkGlobalIds") == 0)
    {
      inScalars = input->GetRowData()->GetGlobalIds();
    }
    else if (strcmp(sel->GetSelectionList()->GetName(), "vtkIndices") == 0)
    {
      use_ids = true;
    }
    else
    {
      inScalars = input->GetRowData()->GetArray(sel->GetSelectionList()->GetName());
    }
  }

  if (inScalars == nullptr && !use_ids)
  {
    vtkErrorMacro("Could not figure out what array to threshold in.");
    return 1;
  }

  int inverse = 0;
  if (sel->GetProperties()->Has(vtkSelectionNode::INVERSE()))
  {
    inverse = sel->GetProperties()->Get(vtkSelectionNode::INVERSE());
  }

  int passThrough = 0;
  if (this->PreserveTopology)
  {
    passThrough = 1;
  }

  int comp_no = 0;
  if (sel->GetProperties()->Has(vtkSelectionNode::COMPONENT_NUMBER()))
  {
    comp_no = sel->GetProperties()->Get(vtkSelectionNode::COMPONENT_NUMBER());
  }

  vtkDataSetAttributes* inRD = input->GetRowData();
  vtkDataSetAttributes* outRD = output->GetRowData();
  vtkSmartPointer<vtkSignedCharArray> rowInArray;
  vtkSmartPointer<vtkIdTypeArray> originalRowIds;
  vtkIdType numRows = input->GetNumberOfRows();

  signed char flag = inverse ? 1 : -1;

  if (passThrough)
  {
    output->ShallowCopy(input);

    rowInArray = vtkSmartPointer<vtkSignedCharArray>::New();
    rowInArray->SetNumberOfComponents(1);
    rowInArray->SetNumberOfTuples(numRows);
    std::fill(rowInArray->GetPointer(0), rowInArray->GetPointer(0) + numRows, flag);
    rowInArray->SetName("vtkInsidedness");
    outRD->AddArray(rowInArray);
  }
  else
  {
    outRD->CopyGlobalIdsOn();
    outRD->CopyAllocate(inRD);

    originalRowIds = vtkSmartPointer<vtkIdTypeArray>::New();
    originalRowIds->SetNumberOfComponents(1);
    originalRowIds->SetName("vtkOriginalRowIds");
    originalRowIds->Allocate(numRows);
    outRD->AddArray(originalRowIds);
  }

  flag = -flag;

  vtkIdType outRCnt = 0;
  for (vtkIdType rowId = 0; rowId < numRows; rowId++)
  {
    int keepRow = this->EvaluateValue(inScalars, comp_no, rowId, lims);
    if (keepRow ^ inverse)
    {
      if (passThrough)
      {
        rowInArray->SetValue(rowId, flag);
      }
      else
      {
        outRD->CopyData(inRD, rowId, outRCnt);
        originalRowIds->InsertNextValue(rowId);
        outRCnt++;
      }
    }
  }
  outRD->Squeeze();
  return 1;
}

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

namespace
{
struct TestItem
{
  template <class ArrayT>
  void operator()(ArrayT* lims, double value, int& keep) const
  {
    auto limsRange = vtk::DataArrayValueRange(lims);
    assert(limsRange.size() % 2 == 0);
    for (auto i = limsRange.begin(); i < limsRange.end(); i += 2)
    {
      if (value >= *i && value <= *(i + 1))
      {
        keep = true;
        return;
      }
    }
    keep = false;
  }

  template <class ArrayT>
  void operator()(
    ArrayT* lims, double value, int& keepCell, int& above, int& below, int& inside) const
  {
    auto limsRange = vtk::DataArrayValueRange(lims);
    assert(limsRange.size() % 2 == 0);
    for (auto i = limsRange.begin(); i < limsRange.end(); i += 2)
    {
      const auto& low = *i;
      const auto& high = *(i + 1);
      if (value >= low && value <= high)
      {
        keepCell = true;
        ++inside;
      }
      else if (value < low)
      {
        ++below;
      }
      else if (value > high)
      {
        ++above;
      }
    }
    keepCell = false;
  }
};
} // namespace

//----------------------------------------------------------------------------
int vtkExtractSelectedThresholds::EvaluateValue(
  vtkDataArray* scalars, int comp_no, vtkIdType id, vtkDataArray* lims)
{
  int keepCell = 0;
  // check the value in the array against all of the thresholds in lims
  // if it is inside any, return true
  double value = 0.0;
  if (comp_no < 0 && scalars)
  {
    // use magnitude.
    int numComps = scalars->GetNumberOfComponents();
    const double* tuple = scalars->GetTuple(id);
    for (int cc = 0; cc < numComps; cc++)
    {
      value += tuple[cc] * tuple[cc];
    }
    value = sqrt(value);
  }
  else
  {
    value = scalars ? scalars->GetComponent(id, comp_no)
                    : static_cast<double>(id); /// <=== precision loss when using id.
  }

  if (!vtkArrayDispatch::Dispatch::Execute(lims, TestItem{}, value, keepCell))
  {
    TestItem{}(lims, value, keepCell);
  }
  return keepCell;
}

//----------------------------------------------------------------------------
int vtkExtractSelectedThresholds::EvaluateValue(vtkDataArray* scalars, int comp_no, vtkIdType id,
  vtkDataArray* lims, int* AboveCount, int* BelowCount, int* InsideCount)
{
  double value = 0.0;
  if (comp_no < 0 && scalars)
  {
    // use magnitude.
    int numComps = scalars->GetNumberOfComponents();
    const double* tuple = scalars->GetTuple(id);
    for (int cc = 0; cc < numComps; cc++)
    {
      value += tuple[cc] * tuple[cc];
    }
    value = sqrt(value);
  }
  else
  {
    value = scalars ? scalars->GetComponent(id, comp_no)
                    : static_cast<double>(id); /// <=== precision loss when using id.
  }

  int keepCell = 0;
  // check the value in the array against all of the thresholds in lims
  // if it is inside any, return true
  int above = 0;
  int below = 0;
  int inside = 0;

  if (!vtkArrayDispatch::Dispatch::Execute(lims, TestItem{}, value, keepCell, above, below, inside))
  {
    TestItem{}(lims, value, keepCell, above, below, inside);
  }

  if (AboveCount)
    *AboveCount = above;
  if (BelowCount)
    *BelowCount = below;
  if (InsideCount)
    *InsideCount = inside;
  return keepCell;
}
