Volume contour spectrum (vtkUnstructuredGrid)

Graphics/CMakeLists.txt

@@ -16,7 +16,7 @@ vtkAppendPolyData.cxx
 vtkAppendSelection.cxx
 vtkApproximatingSubdivisionFilter.cxx
 vtkArcSource.cxx
-vtkAreaContourSignatureFilter.cxx
+vtkAreaContourSpectrumFilter.cxx
 vtkArrayCalculator.cxx
 vtkArrowSource.cxx
 vtkAssignAttribute.cxx
@@ -275,6 +275,7 @@ vtkUnstructuredGridToReebGraphFilter.cxx
 vtkVectorDot.cxx
 vtkVectorNorm.cxx
 vtkVertexGlyphFilter.cxx
+vtkVolumeContourSpectrumFilter.cxx
 vtkVoxelContoursToSurfaceFilter.cxx
 vtkWarpLens.cxx
 vtkWarpScalar.cxx

Graphics/vtkAreaContourSignatureFilter.cxx → Graphics/vtkAreaContourSpectrumFilter.cxx

@@ -12,16 +12,16 @@
      PURPOSE.  See the above copyright notice for more information.
 
 =========================================================================*/
-#include  "vtkAreaContourSignatureFilter.h"
+#include  "vtkAreaContourSpectrumFilter.h"
 
 #include  "vtkInformation.h"
 #include  "vtkInformationVector.h"
 
-vtkCxxRevisionMacro(vtkAreaContourSignatureFilter, "$Revision$");
-vtkStandardNewMacro(vtkAreaContourSignatureFilter);
+vtkCxxRevisionMacro(vtkAreaContourSpectrumFilter, "$Revision$");
+vtkStandardNewMacro(vtkAreaContourSpectrumFilter);
 
 //----------------------------------------------------------------------------
-vtkAreaContourSignatureFilter::vtkAreaContourSignatureFilter()
+vtkAreaContourSpectrumFilter::vtkAreaContourSpectrumFilter()
 {
   this->SetNumberOfInputPorts(2);
   this->ArcId = 0;
@@ -30,12 +30,12 @@ vtkAreaContourSignatureFilter::vtkAreaContourSignatureFilter()
 }
 
 //----------------------------------------------------------------------------
-vtkAreaContourSignatureFilter::~vtkAreaContourSignatureFilter()
+vtkAreaContourSpectrumFilter::~vtkAreaContourSpectrumFilter()
 {
 }
 
 //----------------------------------------------------------------------------
-int vtkAreaContourSignatureFilter::FillInputPortInformation(
+int vtkAreaContourSpectrumFilter::FillInputPortInformation(
   int portNumber, vtkInformation *info)
 {
   switch(portNumber)
@@ -53,7 +53,7 @@ int vtkAreaContourSignatureFilter::FillInputPortInformation(
 }
 
 //----------------------------------------------------------------------------
-int vtkAreaContourSignatureFilter::FillOutputPortInformation(int portNumber,
+int vtkAreaContourSpectrumFilter::FillOutputPortInformation(int portNumber,
   vtkInformation *info)
 {
 
@@ -63,19 +63,19 @@ int vtkAreaContourSignatureFilter::FillOutputPortInformation(int portNumber,
 }
 
 //----------------------------------------------------------------------------
-void vtkAreaContourSignatureFilter::PrintSelf(ostream& os, vtkIndent indent)
+void vtkAreaContourSpectrumFilter::PrintSelf(ostream& os, vtkIndent indent)
 {
   this->Superclass::PrintSelf(os,indent);
 }
 
 //----------------------------------------------------------------------------
-vtkTable* vtkAreaContourSignatureFilter::GetOutput()
+vtkTable* vtkAreaContourSpectrumFilter::GetOutput()
 {
   return vtkTable::SafeDownCast(this->GetOutputDataObject(0));
 }
 
 //----------------------------------------------------------------------------
-int vtkAreaContourSignatureFilter::RequestData(
+int vtkAreaContourSpectrumFilter::RequestData(
   vtkInformation *request, vtkInformationVector **inputVector,
   vtkInformationVector *outputVector)
 {

Graphics/vtkAreaContourSignatureFilter.h → Graphics/vtkAreaContourSpectrumFilter.h

@@ -12,7 +12,7 @@
      PURPOSE.  See the above copyright notice for more information.
 
 =========================================================================*/
-// .NAME vtkAreaContourSignatureFilter - compute an approximation of the area
+// .NAME vtkAreaContourSpectrumFilter - compute an approximation of the area
 // contour signature (evolution of the area of the input surface along an arc of
 // the Reeb graph).
 // .SECTION Description
@@ -33,21 +33,20 @@
 // filter (with customized metrics). It also shows typical vtkReebGraph
 // traversals.
 
-#ifndef __vtkAreaContourSignatureFilter_h
-#define __vtkAreaContourSignatureFilter_h
+#ifndef __vtkAreaContourSpectrumFilter_h
+#define __vtkAreaContourSpectrumFilter_h
 
 #include  "vtkDataObjectAlgorithm.h"
-#include  "vtkDoubleArray.h"
 #include  "vtkTable.h"
 #include  "vtkTriangle.h"
 #include  "vtkReebGraph.h"
 
-class VTK_GRAPHICS_EXPORT vtkAreaContourSignatureFilter :
+class VTK_GRAPHICS_EXPORT vtkAreaContourSpectrumFilter :
   public vtkDataObjectAlgorithm
 {
 public:
-  static vtkAreaContourSignatureFilter* New();
-  vtkTypeRevisionMacro(vtkAreaContourSignatureFilter, vtkDataObjectAlgorithm);
+  static vtkAreaContourSpectrumFilter* New();
+  vtkTypeRevisionMacro(vtkAreaContourSpectrumFilter, vtkDataObjectAlgorithm);
   void PrintSelf(ostream& os, vtkIndent indent);
 
   // Description:
@@ -71,8 +70,8 @@ public:
   vtkTable* GetOutput();
 
 protected:
-  vtkAreaContourSignatureFilter();
-  ~vtkAreaContourSignatureFilter();
+  vtkAreaContourSpectrumFilter();
+  ~vtkAreaContourSpectrumFilter();
 
   vtkIdType ArcId, FieldId;
   int NumberOfSamples;
@@ -84,10 +83,10 @@ protected:
     vtkInformationVector **inputVector, vtkInformationVector *outputVector);
 
 private:
-  vtkAreaContourSignatureFilter(
-    const vtkAreaContourSignatureFilter&);
+  vtkAreaContourSpectrumFilter(
+    const vtkAreaContourSpectrumFilter&);
   // Not implemented.
-  void operator=(const vtkAreaContourSignatureFilter&);
+  void operator=(const vtkAreaContourSpectrumFilter&);
   // Not implemented.
 };
 

Graphics/vtkVolumeContourSpectrumFilter.cxx

+/*=========================================================================
+
+  Program:   Visualization Toolkit
+  Module:    $RCSfile$
+
+  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  "vtkInformation.h"
+#include  "vtkInformationVector.h"
+
+vtkCxxRevisionMacro(vtkVolumeContourSpectrumFilter, "$Revision$");
+vtkStandardNewMacro(vtkVolumeContourSpectrumFilter);
+
+//----------------------------------------------------------------------------
+vtkVolumeContourSpectrumFilter::vtkVolumeContourSpectrumFilter()
+{
+  this->SetNumberOfInputPorts(2);
+  this->ArcId = 0;
+  this->FieldId = 0;
+  this->NumberOfSamples = 100;
+}
+
+//----------------------------------------------------------------------------
+vtkVolumeContourSpectrumFilter::~vtkVolumeContourSpectrumFilter()
+{
+}
+
+//----------------------------------------------------------------------------
+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 portNumber,
+  vtkInformation *info)
+{
+
+  info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkTable");
+
+  return 1;
+}
+
+//----------------------------------------------------------------------------
+void vtkVolumeContourSpectrumFilter::PrintSelf(ostream& os, vtkIndent indent)
+{
+  this->Superclass::PrintSelf(os,indent);
+}
+
+//----------------------------------------------------------------------------
+vtkTable* vtkVolumeContourSpectrumFilter::GetOutput()
+{
+  return vtkTable::SafeDownCast(this->GetOutputDataObject(0));
+}
+
+//----------------------------------------------------------------------------
+int vtkVolumeContourSpectrumFilter::RequestData(
+  vtkInformation *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 = vtkVariantArray::SafeDownCast(
+        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 = vtkDataArray::SafeDownCast(
+        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(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(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);
+      int pos = 0;
+      for(int i = 0; i < NumberOfSamples; i++)
+        {
+        samples[i].first = 0;
+        samples[i].second = 0;
+        while((scalarValues[pos] < min +
+           (i+1.0)*((max - min)/((double)NumberOfSamples)))
+           &&(pos < scalarValues.size()))
+          {
+            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 = i; 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(samples.size());
+      for(int i = 0; i < samples.size(); i++)
+        {
+        outputSignature->SetValue(i, samples[i].second);
+        }
+      output->Initialize();
+      output->AddColumn(outputSignature);
+      outputSignature->Delete();
+      }
+
+    return 1;
+    }
+  return 0;
+}

Graphics/vtkVolumeContourSpectrumFilter.h

+/*=========================================================================
+
+  Program:   Visualization Toolkit
+  Module:    $RCSfile$
+
+  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.
+
+=========================================================================*/
+// .NAME vtkVolumeContourSpectrumFilter - compute an approximation of the
+// volume contour signature (evolution of the volume of the input tet-mesh
+// along an arc of the Reeb graph).
+// .SECTION Description
+// The filter takes a vtkUnstructuredGrid as an input (port 0), along with a
+// vtkReebGraph (port 1).
+// The Reeb graph arc to consider can be specified with SetArcId() (default: 0).
+// The number of (evenly distributed) samples of the signature can be defined
+// with SetNumberOfSamples() (default value: 100).
+// The filter will first try to pull as a scalar field the vtkDataArray with Id
+// 'FieldId' of the vtkUnstructuredGrid, see SetFieldId (default: 0). The
+// filter will abort if this field does not exist.
+//
+// The filter outputs a vtkTable with the volume contour signature
+// approximation, each sample being evenly distributed in the function span of
+// the arc.
+//
+// This filter is a typical example for designing your own contour signature
+// filter (with customized metrics). It also shows typical vtkReebGraph
+// traversals.
+
+#ifndef __vtkVolumeContourSpectrumFilter_h
+#define __vtkVolumeContourSpectrumFilter_h
+
+#include  "vtkDataObjectAlgorithm.h"
+#include  "vtkMath.h"
+#include  "vtkTable.h"
+#include  "vtkTetra.h"
+#include  "vtkReebGraph.h"
+
+class VTK_GRAPHICS_EXPORT vtkVolumeContourSpectrumFilter :
+  public vtkDataObjectAlgorithm
+{
+public:
+  static vtkVolumeContourSpectrumFilter* New();
+  vtkTypeRevisionMacro(vtkVolumeContourSpectrumFilter, vtkDataObjectAlgorithm);
+  void PrintSelf(ostream& os, vtkIndent indent);
+
+  // Description:
+  // Set the arc Id for which the contour signature has to be computed.
+  // Default value: 0
+  vtkSetMacro(ArcId, vtkIdType);
+  vtkGetMacro(ArcId, vtkIdType);
+
+  // Description:
+  // Set the number of samples in the output signature
+  // Default value: 100
+  vtkSetMacro(NumberOfSamples, int);
+  vtkGetMacro(NumberOfSamples, int);
+
+  // Description:
+  // Set the scalar field Id
+  // Default value: 0
+  vtkSetMacro(FieldId, vtkIdType);
+  vtkGetMacro(FieldId, vtkIdType);
+
+  vtkTable* GetOutput();
+
+protected:
+  vtkVolumeContourSpectrumFilter();
+  ~vtkVolumeContourSpectrumFilter();
+
+  vtkIdType ArcId, FieldId;
+  int NumberOfSamples;
+
+  int FillInputPortInformation(int portNumber, vtkInformation *);
+  int FillOutputPortInformation(int portNumber, vtkInformation *info);
+
+  int RequestData(vtkInformation *request,
+    vtkInformationVector **inputVector, vtkInformationVector *outputVector);
+
+private:
+  vtkVolumeContourSpectrumFilter(
+    const vtkVolumeContourSpectrumFilter&);
+  // Not implemented.
+  void operator=(const vtkVolumeContourSpectrumFilter&);
+  // Not implemented.
+};
+
+#endif