Unstructured Gradient Filter and Marching Cubes

[ShaneKirkbride]

I’m trying to show the gradient of a polyData surface similiar to the marching cubes example:
https://lorensen.github.io/VTKExamples/site/Cxx/Modelling/MarchingCubes/

I read in the documentation:
https://vtk.org/doc/nightly/html/classvtkGradientFilter.html
about the unstructured data grid option but I’m unlcear how to use it. Could I get a code snippet or pointer to an example

I saw this but it uses an image:
https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/GradientFilter/

Here is a sample of the code I tried:

#include <string>
#include <vtkSmartPointer.h>
#include <vtkMarchingCubes.h>
#include <vtkVoxelModeller.h>
#include <vtkSphereSource.h>
#include <vtkImageData.h>
#include <vtkDICOMImageReader.h>
#include <vtkJPEGReader.h>
#include <vtkActor.h>
#include <vtkPolyDataMapper.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderWindow.h>
#include <vtkRenderer.h>
#include <vtkImageViewer2.h>

#include <vtkSmartPointer.h>

#include <vtkActor.h>
#include <vtkArrowSource.h>
#include <vtkAssignAttribute.h>
#include <vtkCamera.h>
#include <vtkExtractEdges.h>
#include <vtkGlyph3D.h>
#include <vtkGradientFilter.h>
#include <vtkPolyDataMapper.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkTubeFilter.h>
#include <vtkUnstructuredGridReader.h>

using namespace std;
int main(int argc, char *argv[])
{
  vtkSmartPointer<vtkImageData> volume =
    vtkSmartPointer<vtkImageData>::New();
  double isoValue;
  vtkSmartPointer<vtkSphereSource> sphereSource =
      vtkSmartPointer<vtkSphereSource>::New();
  sphereSource->SetPhiResolution(20);
  sphereSource->SetThetaResolution(20);
  sphereSource->Update();

  double bounds[6];
  sphereSource->GetOutput()->GetBounds(bounds);
  for (unsigned int i = 0; i < 6; i += 2)
  {
      double range = bounds[i + 1] - bounds[i];
      bounds[i] = bounds[i] - .1 * range;
      bounds[i + 1] = bounds[i + 1] + .1 * range;
  }
  vtkSmartPointer<vtkVoxelModeller> voxelModeller =
      vtkSmartPointer<vtkVoxelModeller>::New();
  voxelModeller->SetSampleDimensions(50, 50, 50);
  voxelModeller->SetModelBounds(bounds);
  voxelModeller->SetScalarTypeToFloat();
  voxelModeller->SetMaximumDistance(.1);

  voxelModeller->SetInputConnection(sphereSource->GetOutputPort());
  voxelModeller->Update();
  isoValue = 0.5;
  volume->DeepCopy(voxelModeller->GetOutput());
    
    //std::string inputFilename = argv[1];
    vtkSmartPointer<vtkMarchingCubes> surface =
        vtkSmartPointer<vtkMarchingCubes>::New();
    surface->SetInputData(volume);
    surface->ComputeNormalsOn();
    surface->SetValue(0, isoValue);

    // Create polydata from iso-surface
    vtkSmartPointer<vtkPolyData> marched = vtkSmartPointer<vtkPolyData>::New();
    surface->Update();
    marched->DeepCopy(surface->GetOutput());
    std::cout << "Number of points: " << marched->GetNumberOfPoints() << std::endl;
         
  vtkSmartPointer<vtkExtractEdges> edges = vtkSmartPointer<vtkExtractEdges>::New();
  edges->SetInputData(marched); // ->SetInputData(outputMesh); // ->SetInputConnection(reader->GetOutputPort());

  vtkSmartPointer<vtkTubeFilter> tubes = vtkSmartPointer<vtkTubeFilter>::New();
  tubes->SetInputConnection(edges->GetOutputPort());
  tubes->SetRadius(0.01);
  tubes->SetVaryRadiusToVaryRadiusOff();
  tubes->SetNumberOfSides(4);

  vtkSmartPointer<vtkPolyDataMapper> tubesMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
  tubesMapper->SetInputConnection(tubes->GetOutputPort());
  tubesMapper->SetScalarRange(0.0, 25.0);

  vtkSmartPointer<vtkActor> tubesActor = vtkSmartPointer<vtkActor>::New();
  tubesActor->SetMapper(tubesMapper);

  vtkSmartPointer<vtkGradientFilter> gradients = vtkSmartPointer<vtkGradientFilter>::New();
  gradients->SetInputData(marched); // ->SetInputConnection(reader->GetOutputPort());
  gradients->Update();

  marched->Print(std::cout);
  vtkSmartPointer<vtkAssignAttribute> vectors = vtkSmartPointer<vtkAssignAttribute>::New();
  vectors->SetInputConnection(gradients->GetOutputPort());
  vectors->Assign("Gradients", vtkDataSetAttributes::VECTORS, vtkAssignAttribute::POINT_DATA);

  vtkSmartPointer<vtkArrowSource> arrow = vtkSmartPointer<vtkArrowSource>::New();

  vtkSmartPointer<vtkGlyph3D> glyphs = vtkSmartPointer<vtkGlyph3D>::New();
  glyphs->SetInputConnection(0, vectors->GetOutputPort());
  glyphs->SetInputConnection(1, arrow->GetOutputPort());
  glyphs->ScalingOn();
  glyphs->SetScaleModeToScaleByVector();
  glyphs->SetScaleFactor(0.25);
  glyphs->OrientOn();
  glyphs->ClampingOff();
  glyphs->SetVectorModeToUseVector();
  glyphs->SetIndexModeToOff();

  vtkSmartPointer<vtkPolyDataMapper> glyphMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
  glyphMapper->SetInputConnection(glyphs->GetOutputPort());
  glyphMapper->ScalarVisibilityOff();

  vtkSmartPointer<vtkActor> glyphActor = vtkSmartPointer<vtkActor>::New();
  glyphActor->SetMapper(glyphMapper);

  vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
  renderer->AddActor(tubesActor);
  renderer->AddActor(glyphActor);
  renderer->SetBackground(0.328125, 0.347656, 0.425781);

  vtkSmartPointer<vtkRenderWindow> renwin =
      vtkSmartPointer<vtkRenderWindow>::New();
  renwin->AddRenderer(renderer);
  renwin->SetSize(600, 600);

  renderer->ResetCamera();
  vtkCamera* camera = renderer->GetActiveCamera();
  camera->Elevation(-80.0);
  camera->OrthogonalizeViewUp();
  camera->Azimuth(135.0);

  vtkSmartPointer<vtkRenderWindowInteractor> iren =
      vtkSmartPointer<vtkRenderWindowInteractor>::New();
  iren->SetRenderWindow(renwin);
  renwin->Render();
  iren->Initialize();
  iren->Start();



  return EXIT_SUCCESS;
}

Here is the result - I was hoping for something a little more colorful; and I think I’m probably doing something wrong.
an example like this: Gradient of Unstrucutred Grid in Python
in C++ would probably be helpful enough.

[ShaneKirkbride]

I added some debug values to try to understand this more. I also changed the shape. It seems like all the gradient tuples are 0. I would think they would show some type of change. But I do understand why the shape is all red. I think there is something I’m missing in my pipeline but I’m not sure what. I’ll continue to search. Here is the result and

Here is my code:

#include <string>
#include <vtkActor.h>
#include <vtkArrowSource.h>
#include <vtkAssignAttribute.h>
#include <vtkCamera.h>
#include <vtkConeSource.h>
#include <vtkDoubleArray.h>
#include <vtkSmartPointer.h>
#include <vtkMarchingCubes.h>
#include <vtkVoxelModeller.h>
#include <vtkImageData.h>
#include <vtkImageViewer2.h>
#include <vtkDICOMImageReader.h>
#include <vtkJPEGReader.h>

#include <vtkPolyDataMapper.h>
#include <vtkSphereSource.h>
#include <vtkSmartPointer.h>
#include <vtkExtractEdges.h>
#include <vtkGlyph3D.h>
#include <vtkGradientFilter.h>
#include <vtkPoints.h>
#include <vtkPointData.h>
#include <vtkPolyDataMapper.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkTubeFilter.h>
#include <vtkUnstructuredGridReader.h>

using namespace std;
int main(int argc, char *argv[])
{
    //Create a cone
    vtkSmartPointer<vtkConeSource> coneSource = vtkSmartPointer<vtkConeSource>::New();
    coneSource->Update();
  vtkSmartPointer<vtkImageData> volume = vtkSmartPointer<vtkImageData>::New();
  double isoValue;
  //vtkSmartPointer<vtkSphereSource> sphereSource = vtkSmartPointer<vtkSphereSource>::New();
  //sphereSource->SetPhiResolution(200);
  //sphereSource->SetThetaResolution(200);
  //sphereSource->Update();

  double bounds[6];
  coneSource->GetOutput()->GetBounds(bounds);
  for (unsigned int i = 0; i < 6; i += 2)
  {
      double range = bounds[i + 1] - bounds[i];
      bounds[i] = bounds[i] - .1 * range;
      bounds[i + 1] = bounds[i + 1] + .1 * range;
  }
  vtkSmartPointer<vtkVoxelModeller> voxelModeller =
      vtkSmartPointer<vtkVoxelModeller>::New();
  voxelModeller->SetSampleDimensions(50, 50, 50);
  voxelModeller->SetModelBounds(bounds);
  voxelModeller->SetScalarTypeToFloat();
  voxelModeller->SetMaximumDistance(.1);

  voxelModeller->SetInputConnection(coneSource->GetOutputPort());
  voxelModeller->Update();
  isoValue = 0.5;
  volume->DeepCopy(voxelModeller->GetOutput());
    
    //std::string inputFilename = argv[1];
    vtkSmartPointer<vtkMarchingCubes> surface =
        vtkSmartPointer<vtkMarchingCubes>::New();
    surface->SetInputData(volume);
    surface->ComputeNormalsOn();
    surface->SetValue(0, isoValue);

    // Create polydata from iso-surface
    vtkSmartPointer<vtkPolyData> marched = vtkSmartPointer<vtkPolyData>::New();
    surface->Update();
    marched->DeepCopy(surface->GetOutput());
    std::cout << "Number of points: " << marched->GetNumberOfPoints() << std::endl;
 
  vtkSmartPointer<vtkGradientFilter> gradientFilter = vtkSmartPointer<vtkGradientFilter>::New();
  gradientFilter->SetInputData(marched); //->SetInputConnection(volume->GetOutputPort());
  gradientFilter->Update();
  
  auto gradientOutput = gradientFilter->GetOutput();
  gradientOutput->Print(std::cout);
  ////save gradient array
  vtkSmartPointer<vtkDoubleArray> gradients = vtkSmartPointer<vtkDoubleArray>::New();

  gradients->SetNumberOfTuples(gradientOutput->GetNumberOfPoints());
  gradients->SetName("gradients");
  gradients->Print(std::cout); 
  
  int numberOfTuples = gradientOutput->GetNumberOfPoints();
  auto gradient_array = gradientOutput->GetPointData()->GetArray("Gradients");
  for (vtkIdType ii = 0; ii < numberOfTuples; ii++) {
      double* gradientValue = gradient_array->GetTuple(ii);

      std::cout << gradientValue[0] << " " 
                << gradientValue[1] << " "
                << gradientValue[2] << std::endl;

  }
  
  vtkSmartPointer<vtkExtractEdges> edges = vtkSmartPointer<vtkExtractEdges>::New();
  edges->SetInputData(gradientOutput); // ->SetInputData(outputMesh); // ->SetInputConnection(reader->GetOutputPort());

  vtkSmartPointer<vtkTubeFilter> tubes = vtkSmartPointer<vtkTubeFilter>::New();
  tubes->SetInputConnection(edges->GetOutputPort());
  tubes->SetRadius(0.01);
  tubes->SetVaryRadiusToVaryRadiusOff();
  tubes->SetNumberOfSides(4);

  vtkSmartPointer<vtkPolyDataMapper> tubesMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
  tubesMapper->SetInputConnection(tubes->GetOutputPort());
  tubesMapper->SetScalarRange(0.0, 25.0);

  vtkSmartPointer<vtkActor> tubesActor = vtkSmartPointer<vtkActor>::New();
  tubesActor->SetMapper(tubesMapper);
  
  vtkSmartPointer<vtkAssignAttribute> vectors = vtkSmartPointer<vtkAssignAttribute>::New();
  vectors->SetInputConnection(gradientFilter->GetOutputPort());
  vectors->Assign("Gradients", vtkDataSetAttributes::VECTORS, vtkAssignAttribute::POINT_DATA);

  vtkSmartPointer<vtkArrowSource> arrow = vtkSmartPointer<vtkArrowSource>::New();

  vtkSmartPointer<vtkGlyph3D> glyphs = vtkSmartPointer<vtkGlyph3D>::New();
  glyphs->SetInputConnection(0, vectors->GetOutputPort());
  glyphs->SetInputConnection(1, arrow->GetOutputPort());
  glyphs->ScalingOn();
  glyphs->SetScaleModeToScaleByVector();
  glyphs->SetScaleFactor(0.25);
  glyphs->OrientOn();
  glyphs->ClampingOff();
  glyphs->SetVectorModeToUseVector();
  glyphs->SetIndexModeToOff();

  vtkSmartPointer<vtkPolyDataMapper> glyphMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
  glyphMapper->SetInputConnection(glyphs->GetOutputPort());
  glyphMapper->ScalarVisibilityOff();

  vtkSmartPointer<vtkActor> glyphActor = vtkSmartPointer<vtkActor>::New();
  glyphActor->SetMapper(glyphMapper);

  vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
  renderer->AddActor(tubesActor);
  renderer->AddActor(glyphActor);
  renderer->SetBackground(0.328125, 0.347656, 0.425781);

  vtkSmartPointer<vtkRenderWindow> renwin =
      vtkSmartPointer<vtkRenderWindow>::New();
  renwin->AddRenderer(renderer);
  renwin->SetSize(600, 600);

  renderer->ResetCamera();
  vtkCamera* camera = renderer->GetActiveCamera();
  camera->Elevation(-80.0);
  camera->OrthogonalizeViewUp();
  camera->Azimuth(135.0);

  vtkSmartPointer<vtkRenderWindowInteractor> iren =
      vtkSmartPointer<vtkRenderWindowInteractor>::New();
  iren->SetRenderWindow(renwin);
  renwin->Render();
  iren->Initialize();
  iren->Start();



  return EXIT_SUCCESS;
}

[ShaneKirkbride]

Can someone please help me use the gradient filter on an isoSurface?