Writing custom class based on `vtkDataSet`

Hi,

I’m thinking to write custom volume dataset for working with scalars.
I work with geo data and usually there we have arbitrary X,Y coordinates (XY are oriented in horizontal plane, Z is vertical axis) and Z is usually have constant length and constant spacing.

For example for for some points X0, Y0, we have Z array [10, 8, 6, 4, 2, 0, -2,-4, -6, -8].
For X1, Y1 we may have Z array [2, 0, -2, -4, -6, -8, -10, -12, -14, -16]. Here is the schema of that:
image

At both surface coordinates (X0Y0 and X1Y1) Z have constant spacing (-2) and constant length (10 elements). In the same time X and Y coordinates are absolutely arbitrary.

I plan to write a class similar to vtkImageData but adjusted to my special case, thus reducing memory usage. Because my data is pretty big (may take few gigabytes) and I don’t want to store every XYZ coordinate for each point, they can be calculated pretty simple.

The visualization should be similar to vtkImageData based on scalars.

As I don’t have an experience of writing vtkDataSet classes I need an advice. Maybe the idea is completely absurd?

No this is not absurd depending on the performance you want to wring out of it, but it is not for the faint of heart :slight_smile: You also have to consider long-term maintenance, which is a burden that falls on the entire VTK community, so it’s important that there is enough value to warrant the addition. Of course maybe you are thinking of creating a separate repo, with appropriate CMake files, that is meant to work in conjunction with VTK (that you/your collaborators maintain and test etc).

If you decide to create a new dataset, many of the vtkDataSet methods are not too bad to implement. Since you have a hybrid implicit + explicit representation, you have to figure out how to number the points and cells so that GetPoint() and GetCell() make sense, and how to create (on the fly) concrete cells and points that are implicitly represented in your dataset. FindCell() and FindPoint() are probably some of the trickiest methods.

It wouldn’t be too hard initially to prototype using a reader that transforms a file of xy-points + 10-component scalars into an unstructured grid. (I assume that you have some sort of surface tessellation that goes along with the x-y points? And that these tessellated cells are triangles or quads that when extruded create 3D wedges or hexes?)

Have fun!

1 Like

Thank you for the response,

Absolutely right. About an hour ago I finally achieved some success (with memory leaks yet of course) and I’m proud to show it :slight_smile:

image

As you said one have to provide XYZ surface grid and Z spacings (double array where all elements are either positive or negative). To visualize with colors o course one need to provide array of values for each cell/point.

Does these methods should return the nearest point/cell in the terms of Eigen space like min(sqrt(x1-x0)^2+(y1-y0)^2+(z1-z0)^2) ? For now I have defined these methods with that logic.

I would like to contribute to VTK but not soon as I need to be confident that the project is pretty stable and may be applied in practical tasks.

Does these methods should return the nearest point/cell in the terms of Eigen space like min(sqrt(x1-x0)^2+(y1-y0)^2+(z1-z0)^2) ? For now I have defined these methods with that logic.

Yes, what is the closest point to a specified position x[3]; and what cell contains x (or is closest to x within a tolerance).

1 Like

@will.schroeder encountered a problem with points indexing and data coloring.

In the picture above you can see that scalar data continuously increases in XY plane (color biggest gradient is oriented along Z axis).

I need to set scalars that are oriented along Z axis and thus the color biggest gradient should be set along X or Y direction. The reason is that usually my data is stored in that way and as it is pretty big (few gigabytes) I don’t want make any manipulation with data (transposition or anything).

Is there a way I could map scalar for points id? For example I guess now for point[0] the color is defined by scalar[0] and for point[9] by scalar[9]. I need to change it for example in the following way: point[9] should be defined by scalar[49]

This is a major design choice of creating a new dataset, you have to define the ordering of points and cells in such a way as it is consistent with the data representation and attributes. It sounds like you need to rethink the point ordering to be consistent with the point attribute ordering (i.e., scalars); there’s no reason that the point ordering cannot run “vertically” from the base plane through the points associated with the scalars, and so on for each point.

1 Like

Thank you a lot,

I needed to correctly set id for each point location. Maybe for someone it will be useful:

  // set points in vertical order
  vtkNew<vtkPoints> points;
  points->InsertNextPoint(0, 0, 0);
  points->InsertNextPoint(0, 0, 1);
  points->InsertNextPoint(1, 0, 0);
  points->InsertNextPoint(1, 0, 1);
  points->InsertNextPoint(0, 1, 0);
  points->InsertNextPoint(0, 1, 1);

  // set the correct id for each point location
  vtkNew<vtkWedge> wedge;
  wedge->GetPointIds()->SetId(0,0);
  wedge->GetPointIds()->SetId(1,2);
  wedge->GetPointIds()->SetId(2,4);
  wedge->GetPointIds()->SetId(3,1);
  wedge->GetPointIds()->SetId(4,3);
  wedge->GetPointIds()->SetId(5,5);

Example that eliminates some uncertainties about scalar-points connection. I leave here in case I forget something:

#include <vtkActor.h>
#include <vtkCamera.h>
#include <vtkCellData.h>
#include <vtkColor.h>
#include <vtkDataSetMapper.h>
#include <vtkDoubleArray.h>
#include <vtkInteractorStyleSwitch.h>
#include <vtkLookupTable.h>
#include <vtkMapper.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkPointData.h>
#include <vtkPoints.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkAxesActor.h>
#include <vtkWedge.h>
#include <vtkTextActor.h>
#include <vtkScalarBarActor.h>
#include <vtkLabeledDataMapper.h>
#include <vtkTextProperty.h>

int main(int, char**) {
  vtkNew<vtkNamedColors> colors;
  // Set the background color.
  std::array<unsigned char, 4> bkg{{26, 51, 102, 255}};
  colors->SetColor("BkgColor", bkg.data());

  int numberOfVertices = 6;


  // Setup wedge points in usual order (one triangle face after another triangle face)
  vtkNew<vtkPoints> normalWedgePoints;
  normalWedgePoints->InsertNextPoint(0, 0, 0);
  normalWedgePoints->InsertNextPoint(1, 0, 0);
  normalWedgePoints->InsertNextPoint(0, 1, 0);
  normalWedgePoints->InsertNextPoint(0, 0, 1);
  normalWedgePoints->InsertNextPoint(1, 0, 1);
  normalWedgePoints->InsertNextPoint(0, 1, 1);

  // Eliminate uncertainty how scalaras are linked with points
  vtkNew<vtkDoubleArray> normalWedgeScalars;
  normalWedgeScalars->SetName("NormalWedgeArray");
  normalWedgeScalars->SetNumberOfComponents(1);
  normalWedgeScalars->SetNumberOfValues(numberOfVertices);
  for (int i = 0; i < numberOfVertices; i++)
  {
    normalWedgeScalars->SetValue(i, i);
  }

  vtkNew<vtkWedge> normalWedge;
  for (int i = 0; i < numberOfVertices; i++)
  {
    normalWedge->GetPointIds()->SetId(i, i);
  }

  vtkNew<vtkUnstructuredGrid> normalWedgeUGrid;
  normalWedgeUGrid->SetPoints(normalWedgePoints);
  normalWedgeUGrid->InsertNextCell(normalWedge->GetCellType(), normalWedge->GetPointIds());
  normalWedgeUGrid->GetPointData()->SetScalars(normalWedgeScalars);


  // Setup wedge points in vertical order
  vtkNew<vtkPoints> verticalWedgePoints;
  verticalWedgePoints->InsertNextPoint(0, 0, 0);
  verticalWedgePoints->InsertNextPoint(0, 0, 1);
  verticalWedgePoints->InsertNextPoint(1, 0, 0);
  verticalWedgePoints->InsertNextPoint(1, 0, 1);
  verticalWedgePoints->InsertNextPoint(0, 1, 0);
  verticalWedgePoints->InsertNextPoint(0, 1, 1);

  // Eliminate uncertainty how scalaras are linked with points
  vtkNew<vtkDoubleArray> verticalWedgeScalars;
  verticalWedgeScalars->SetName("VerticalWedgeArray");
  verticalWedgeScalars->SetNumberOfComponents(1);
  verticalWedgeScalars->SetNumberOfValues(numberOfVertices);
  for (int i = 0; i < numberOfVertices; i++)
  {
    verticalWedgeScalars->SetValue(i, i);
  }

  vtkNew<vtkWedge> verticalWedge;
  verticalWedge->GetPointIds()->SetId(0,0);
  verticalWedge->GetPointIds()->SetId(1,2);
  verticalWedge->GetPointIds()->SetId(2,4);
  verticalWedge->GetPointIds()->SetId(3,1);
  verticalWedge->GetPointIds()->SetId(4,3);
  verticalWedge->GetPointIds()->SetId(5,5);

  vtkNew<vtkUnstructuredGrid> verticalWedgeUGrid;
  verticalWedgeUGrid->SetPoints(verticalWedgePoints);
  verticalWedgeUGrid->InsertNextCell(verticalWedge->GetCellType(), verticalWedge->GetPointIds());
  verticalWedgeUGrid->GetPointData()->SetScalars(verticalWedgeScalars);

  // Setup mappers and actors
  vtkNew<vtkLookupTable> lut;
  lut->SetNumberOfTableValues(6);
  lut->Build();

  vtkNew<vtkDataSetMapper> normalWedgeMapper;
  normalWedgeMapper->SetInputData(normalWedgeUGrid);
  normalWedgeMapper->SetLookupTable(lut);
  normalWedgeMapper->SetScalarRange(0, numberOfVertices-1);
  normalWedgeMapper->ScalarVisibilityOn();
  normalWedgeMapper->SetScalarModeToUsePointData();

  vtkNew<vtkDataSetMapper> verticalWedgeMapper;
  verticalWedgeMapper->SetInputData(verticalWedgeUGrid);
  verticalWedgeMapper->SetLookupTable(lut);
  verticalWedgeMapper->SetScalarRange(0, numberOfVertices-1);
  verticalWedgeMapper->ScalarVisibilityOn();
  verticalWedgeMapper->SetScalarModeToUsePointData();

  vtkNew<vtkActor> normalWedgeActor;
  normalWedgeActor->SetMapper(normalWedgeMapper);
  normalWedgeActor->GetProperty()->EdgeVisibilityOn();
  normalWedgeActor->GetProperty()->SetOpacity(1);

  vtkNew<vtkActor> verticalWedgeActor;
  verticalWedgeActor->SetMapper(verticalWedgeMapper);
  verticalWedgeActor->GetProperty()->EdgeVisibilityOn();
  verticalWedgeActor->GetProperty()->SetOpacity(1);

  vtkNew<vtkScalarBarActor> scalarBar;
  scalarBar->SetLookupTable(lut);
  scalarBar->SetTitle("Values");
  scalarBar->SetNumberOfLabels(6);
  scalarBar->UnconstrainedFontSizeOn();

  vtkNew<vtkTextActor> normalWedgeTextActor;
  normalWedgeTextActor->SetInput("Normal wedge points order");
  normalWedgeTextActor->SetPosition(20, 20);

  vtkNew<vtkTextActor> verticalWedgeTextActor;
  verticalWedgeTextActor->SetInput("Vertical wedge points order");
  verticalWedgeTextActor->SetPosition(20, 20);

  vtkNew<vtkLabeledDataMapper> normalWedgeLabelMapper;
  normalWedgeLabelMapper->SetInputData(normalWedgeUGrid);

  vtkNew<vtkLabeledDataMapper> verticalWedgeLabelMapper;
  verticalWedgeLabelMapper->SetInputData(verticalWedgeUGrid);

  vtkNew<vtkActor2D> normalWedgeLabelActor;
  normalWedgeLabelActor->SetMapper(normalWedgeLabelMapper);

  vtkNew<vtkActor2D> verticalWedgeLabelActor;
  verticalWedgeLabelActor->SetMapper(verticalWedgeLabelMapper);

  vtkNew<vtkAxesActor> axes;

  vtkNew<vtkRenderer> rendererLeft;
  rendererLeft->SetBackground(colors->GetColor3d("BkgColor").GetData());
  rendererLeft->AddActor(normalWedgeActor);
  rendererLeft->AddActor(normalWedgeLabelActor);
  rendererLeft->AddActor(normalWedgeTextActor);
  rendererLeft->AddActor(scalarBar);
  rendererLeft->AddActor(axes);
  rendererLeft->ResetCamera();
  rendererLeft->GetActiveCamera()->Azimuth(210);
  rendererLeft->GetActiveCamera()->Elevation(-30);
  rendererLeft->ResetCameraClippingRange();
  rendererLeft->SetViewport(0, 0, 0.5, 1);

  vtkNew<vtkRenderer> rendererRight;
  rendererRight->SetBackground(colors->GetColor3d("BkgColor").GetData());
  rendererRight->AddActor(verticalWedgeActor);
  rendererRight->AddActor(verticalWedgeLabelActor);
  rendererRight->AddActor(verticalWedgeTextActor);
  rendererRight->AddActor(scalarBar);
  rendererRight->AddActor(axes);
  rendererRight->ResetCamera();
  rendererRight->GetActiveCamera()->Azimuth(210);
  rendererRight->GetActiveCamera()->Elevation(-30);
  rendererRight->ResetCameraClippingRange();
  rendererRight->SetViewport(0.5, 0, 1, 1);

  vtkNew<vtkRenderWindow> renderWindow;
  renderWindow->AddRenderer(rendererLeft);
  renderWindow->AddRenderer(rendererRight);
  renderWindow->SetSize(640, 480);

  vtkNew<vtkRenderWindowInteractor> interactor;
  vtkNew<vtkInteractorStyleSwitch> style;
  interactor->SetInteractorStyle(style);
  interactor->SetRenderWindow(renderWindow);

  renderWindow->SetWindowName("Wedge points ordered normally and vertically");
  renderWindow->Render();

  interactor->Start();
}