/Filters/General/vtkVolumeContourSpectrumFilter.cxx

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

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

#include "vtkDataArray.h"
#include "vtkDataSetAttributes.h"
#include "vtkEdgeListIterator.h"
#include "vtkIdList.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkReebGraph.h"
#include "vtkTable.h"
#include "vtkTetra.h"
#include "vtkUnstructuredGrid.h"
#include "vtkVariantArray.h"

vtkStandardNewMacro(vtkVolumeContourSpectrumFilter);

//------------------------------------------------------------------------------
vtkVolumeContourSpectrumFilter::vtkVolumeContourSpectrumFilter()
{
  this->SetNumberOfInputPorts(2);
  this->ArcId = 0;
  this->FieldId = 0;
  this->NumberOfSamples = 100;
}

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

//------------------------------------------------------------------------------
int vtkVolumeContourSpectrumFilter::FillInputPortInformation(int portNumber, vtkInformation* info)
{
  switch (portNumber)
  {
    case 0:
      info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE());
      info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkUnstructuredGrid");
      break;
    case 1:
      info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE());
      info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkReebGraph");
      break;
  }
  return 1;
}

//------------------------------------------------------------------------------
int vtkVolumeContourSpectrumFilter::FillOutputPortInformation(
  int vtkNotUsed(portNumber), vtkInformation* info)
{

  info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkTable");

  return 1;
}

//------------------------------------------------------------------------------
void vtkVolumeContourSpectrumFilter::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
  os << indent << "Arc Id: " << this->ArcId << "\n";
  os << indent << "Number of Samples: " << this->NumberOfSamples << "\n";
  os << indent << "Field Id: " << this->FieldId << "\n";
}

//------------------------------------------------------------------------------
vtkTable* vtkVolumeContourSpectrumFilter::GetOutput()
{
  return vtkTable::SafeDownCast(this->GetOutputDataObject(0));
}

//------------------------------------------------------------------------------
int vtkVolumeContourSpectrumFilter::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{

  vtkInformation *inInfoMesh = inputVector[0]->GetInformationObject(0),
                 *inInfoGraph = inputVector[1]->GetInformationObject(0);

  if ((!inInfoMesh) || (!inInfoGraph))
  {
    return 0;
  }
  vtkUnstructuredGrid* inputMesh =
    vtkUnstructuredGrid::SafeDownCast(inInfoMesh->Get(vtkUnstructuredGrid::DATA_OBJECT()));

  vtkReebGraph* inputGraph =
    vtkReebGraph::SafeDownCast(inInfoGraph->Get(vtkReebGraph::DATA_OBJECT()));

  if ((inputMesh) && (inputGraph))
  {
    vtkInformation* outInfo = outputVector->GetInformationObject(0);
    vtkTable* output = vtkTable::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

    if (output)
    {

      // Retrieve the arc given by ArcId
      vtkVariantArray* edgeInfo = vtkArrayDownCast<vtkVariantArray>(
        inputGraph->GetEdgeData()->GetAbstractArray("Vertex Ids"));
      // Invalid Reeb graph (no information associated to the edges)
      if (!edgeInfo)
        return 0;

      // Retrieve the information to get the critical vertices Ids
      vtkDataArray* criticalPointIds =
        vtkArrayDownCast<vtkDataArray>(inputGraph->GetVertexData()->GetAbstractArray("Vertex Ids"));
      // Invalid Reeb graph (no information associated to the vertices)
      if (!criticalPointIds)
        return 0;

      vtkAbstractArray* vertexList = edgeInfo->GetPointer(ArcId)->ToArray();

      // the arc defined by ArcId does not exist (out of bound?)
      if (!vertexList)
        return 0;

      vtkDataArray* scalarField = inputMesh->GetPointData()->GetArray(FieldId);
      if (!scalarField)
        return 0;

      // parse the input vertex list (region in which the connectivity of the
      // level sets does not change) and compute the area signature
      double cumulativeVolume = 0;
      std::vector<double> scalarValues, volumeSignature;
      std::vector<int> vertexIds;
      std::vector<bool> visitedTetrahedra;

      visitedTetrahedra.resize(inputMesh->GetNumberOfCells());
      vertexIds.resize(vertexList->GetNumberOfTuples() + 2);
      scalarValues.resize(vertexIds.size());
      volumeSignature.resize(vertexIds.size());

      // include the critical points in the computation
      //  - iterates through the edges of the Reeb graph until we found the arc
      //  we're looking for
      //  - retrieve the Source and Target of the edge
      //  - pick the corresponding mesh vertex Ids in the VertexData.
      std::pair<int, int> criticalPoints;
      vtkEdgeListIterator* eIt = vtkEdgeListIterator::New();
      inputGraph->GetEdges(eIt);

      do
      {
        vtkEdgeType e = eIt->Next();
        if (e.Id == ArcId)
        {
          if ((criticalPointIds->GetTuple(e.Source)) && (criticalPointIds->GetTuple(e.Target)))
          {
            criticalPoints.first = (int)*(criticalPointIds->GetTuple(e.Source));
            criticalPoints.second = (int)*(criticalPointIds->GetTuple(e.Target));
          }
          else
            // invalid Reeb graph
            return 0;
        }
      } while (eIt->HasNext());

      eIt->Delete();

      vertexIds[0] = criticalPoints.first;
      vertexIds[vertexIds.size() - 1] = criticalPoints.second;
      // NB: the vertices of vertexList are already in sorted order of function
      // value.
      for (int i = 0; i < vertexList->GetNumberOfTuples(); i++)
        vertexIds[i + 1] = vertexList->GetVariantValue(i).ToInt();

      // mark all the input triangles as non visited.
      for (unsigned int i = 0; i < visitedTetrahedra.size(); i++)
        visitedTetrahedra[i] = false;

      // now do the parsing
      double min = scalarField->GetComponent(vertexIds[0], 0),
             max = scalarField->GetComponent(vertexIds[vertexIds.size() - 1], 0);
      for (unsigned int i = 0; i < vertexIds.size(); i++)
      {
        scalarValues[i] = scalarField->GetComponent(vertexIds[i], 0);

        vtkIdList* starTetrahedronList = vtkIdList::New();
        inputMesh->GetPointCells(vertexIds[i], starTetrahedronList);

        for (int j = 0; j < starTetrahedronList->GetNumberOfIds(); j++)
        {
          vtkIdType tId = starTetrahedronList->GetId(j);
          if (!visitedTetrahedra[tId])
          {
            vtkTetra* t = vtkTetra::SafeDownCast(inputMesh->GetCell(tId));

            if ((scalarField->GetComponent(t->GetPointIds()->GetId(0), 0) <= scalarValues[i]) &&
              (scalarField->GetComponent(t->GetPointIds()->GetId(1), 0) <= scalarValues[i]) &&
              (scalarField->GetComponent(t->GetPointIds()->GetId(2), 0) <= scalarValues[i]) &&
              (scalarField->GetComponent(t->GetPointIds()->GetId(3), 0) <= scalarValues[i]) &&
              (scalarField->GetComponent(t->GetPointIds()->GetId(0), 0) >= min) &&
              (scalarField->GetComponent(t->GetPointIds()->GetId(1), 0) >= min) &&
              (scalarField->GetComponent(t->GetPointIds()->GetId(2), 0) >= min) &&
              (scalarField->GetComponent(t->GetPointIds()->GetId(3), 0) >= min))
            {
              // make sure the triangle is strictly in the covered function
              // span.

              double p0[3], p1[3], p2[3], p3[3];
              inputMesh->GetPoint(t->GetPointIds()->GetId(0), p0);
              inputMesh->GetPoint(t->GetPointIds()->GetId(1), p1);
              inputMesh->GetPoint(t->GetPointIds()->GetId(2), p2);
              inputMesh->GetPoint(t->GetPointIds()->GetId(3), p3);
              cumulativeVolume += vtkTetra::ComputeVolume(p0, p1, p2, p3);

              visitedTetrahedra[tId] = true;
            }
          }
        }
        volumeSignature[i] = cumulativeVolume;
        starTetrahedronList->Delete();
      }

      // now adjust to the desired sampling
      std::vector<std::pair<int, double>> samples(NumberOfSamples);
      unsigned int pos = 0;
      for (int i = 0; i < NumberOfSamples; i++)
      {
        samples[i].first = 0;
        samples[i].second = 0;
        double temp = min + (i + 1.0) * ((max - min) / ((double)NumberOfSamples));
        while ((pos < scalarValues.size()) && (scalarValues[pos] < temp))
        {
          samples[i].first++;
          samples[i].second += volumeSignature[pos];
          pos++;
        }
        if (samples[i].first)
        {
          samples[i].second /= samples[i].first;
        }
      }

      // no value at the start? put 0
      if (!samples[0].first)
      {
        samples[0].first = 1;
        samples[0].second = 0;
      }
      // no value at the end? put the cumulative area
      if (!samples[samples.size() - 1].first)
      {
        samples[samples.size() - 1].first = 1;
        samples[samples.size() - 1].second = cumulativeVolume;
      }

      // fill out the blanks
      int lastSample = 0;
      for (int i = 0; i < NumberOfSamples; i++)
      {
        if (!samples[i].first)
        {
          // not enough vertices in the region for the number of desired
          // samples. we have to interpolate.

          // first, search for the next valid sample
          int nextSample = i;
          for (; nextSample < NumberOfSamples; nextSample++)
          {
            if (samples[nextSample].first)
              break;
          }

          // next interpolate
          samples[i].second = samples[lastSample].second +
            (i - lastSample) * (samples[nextSample].second - samples[lastSample].second) /
              (nextSample - lastSample);
        }
        else
          lastSample = i;
      }

      // now prepare the output
      vtkVariantArray* outputSignature = vtkVariantArray::New();
      outputSignature->SetNumberOfTuples(static_cast<vtkIdType>(samples.size()));
      for (unsigned int i = 0; i < samples.size(); i++)
      {
        outputSignature->SetValue(i, samples[i].second);
      }
      output->Initialize();
      output->AddColumn(outputSignature);
      outputSignature->Delete();
    }

    return 1;
  }
  return 0;
}