Cell Set as Unstructured Grid

In the last 2 blogs (1, 2), we have looked at creating vtkMappedDataArray<> subclasses to represent different memory layouts. This enabled us to represent various types of grids without explicitly storing the point coordinates. In this blog, we will take this idea to the next step and create an unstructured grid that uses a vtkMappedDataArray<> to represent its points but also store its connectivity implicitly. There are many use cases for this. For example, we can develop a better threshold filter. The current implementation (vtkThreshold) creates an unstructured grid no matter the input type. For structured data, this can be terribly memory inefficient because it creates connectivity arrays and potentially a new explicit point coordinates array (for vtkImageData and vtkRectilinearGrid). One alternative to this approach is to create a new unstructured grid type that stores the input dataset and a set of cell ids that meet the threshold criteria. We will cover how such a class can be defined using the zero copy infrastructure. I will describe the new threshold algorithm in detail in a future blog.

vtkImagePointsArray

First, let's recap the points array from a previous blog. The full source for this array is here and here. In summary, vtkImagePointsArray is a sub-class of vtkMappedArray<> that obtains its values from a vtkImageData rather than from an array stored explicitly. This enables us to represent unstructured grids with implicit point coordinates. The meat of this class is as follows.

template <class Scalar> Scalar& vtkImagePointsArray<Scalar>
::GetValueReference(vtkIdType idx)
{
  assert(this->Points);

  const vtkIdType tuple = idx / 3;
  const vtkIdType comp = idx % 3;

  return this->Points->GetPoint(tuple)[comp];
}

Note that Points is actually a vtkImageData, which calculates point coordinates implicitly with the following code.

void vtkImageData::GetPoint(vtkIdType ptId, double x[3])
{
  int i, loc[3];
  const double *origin = this->Origin;
  const double *spacing = this->Spacing;
  const int* extent = this->Extent;

  vtkIdType dims[3];
  dims[0] = extent[1] - extent[0] + 1;
  dims[1] = extent[3] - extent[2] + 1;
  dims[2] = extent[5] - extent[4] + 1;

  loc[0] = ptId % dims[0];
  loc[1] = (ptId / dims[0]) % dims[1];
  loc[2] = ptId / (dims[0]*dims[1]);

  for (i=0; i<3; i++)
    {
    x[i] = origin[i] + (loc[i]+extent[i*2]) * spacing[i];
    }
}

I left out some of the code for simplicity. We can exercise this class with the following.

vtkNew<vtkRTAnalyticSource> source;
source->SetWholeExtent(-50, 50, -50, 50, -50, 50);
source->Update();

vtkImageData* wavelet = source->GetOutput();

vtkNew<vtkImageData> img;
img->CopyStructure(wavelet);

vtkNew<vtkImagePointsArray<double> > testPts;
testPts->InitializeArray(img.GetPointer());
testPts->SetName("pts");

vtkNew<vtkPoints> points;
points->SetData(testPts.GetPointer());

// This creates an unstructured grid from the input
// image data.
vtkNew<vtkCleanUnstructuredGrid> toUGrid;
toUGrid->SetInputData(wavelet);
toUGrid->Update();

vtkUnstructuredGrid* ugrid = toUGrid->GetOutput();
// Here we substitute the implicit points.
ugrid->SetPoints(points.GetPointer());

vtkNew<vtkContourFilter> contour2;
contour2->SetInputData(ugrid);
contour2->SetValue(0, 200);
contour2->Update();

vtkCellSet

Our next step is to define a new unstructured grid that does not store a new connectivity or point coordinate array for an input grid. Instead, we want it to simply store the ids of the cells selected by the threshold process. We will do this by leveraging the mapped unstructured grid functionality created by David Lonie and described here. We will create a new unstructured grid (subclass of vtkMappedUnstructuredGrid) called vtkCellSet. Our first step is to create a helper class (subclass of vtkObject) that takes care of the implementation of core functions. See here and here for the full implementation.

As I mentioned, this class should store a dataset and a set of selected cells. This looks as follows.

class vtkCellSetImpl : public vtkObject
{
public:
  void SetDataSet(vtkDataSet*);
  void SetCellIds(vtkIdTypeArray*);

private:
  vtkDataSet* DataSet;
  vtkIdTypeArray* CellIds;
};

vtkCxxSetObjectMacro(vtkCellSetImpl,DataSet,vtkDataSet);
vtkCxxSetObjectMacro(vtkCellSetImpl,CellIds,vtkIdTypeArray);

Then we need to implement a few key methods as follows.

vtkIdType vtkCellSetImpl::GetNumberOfCells()
{
  return this->CellIds->GetNumberOfTuples();
}

int vtkCellSetImpl::GetCellType(vtkIdType id)
{
  return this->DataSet->GetCellType(this->CellIds->GetValue(id));
}

void vtkCellSetImpl::GetCellPoints(vtkIdType cellId,
                                   vtkIdList *ptIds)
{
  return this->DataSet->GetCellPoints(this->CellIds->GetValue(cellId),
    ptIds);
}

int vtkCellSetImpl::GetMaxCellSize()
{
  return this->DataSet->GetMaxCellSize();
}

Note that many of these are deferred to the underlying dataset (DataSet). The key point is that cell ids are transformed using the internal selection array (CellIds) as in this example.

void vtkCellSetImpl::GetCellPoints(vtkIdType cellId,
                                   vtkIdList *ptIds)
{
  return this->DataSet->GetCellPoints(this->CellIds->GetValue(cellId),
    ptIds);
}

This means that any DataSet cell not listed in the CellIds array will be ignored by this class. We can exercise this class as follows

vtkNew<vtkRTAnalyticSource> source;
source->SetWholeExtent(-50, 50, -50, 50, -50, 50);
source->Update();

vtkImageData* wavelet = source->GetOutput();

vtkNew<vtkImageData> img;
img->CopyStructure(wavelet);

vtkNew<vtkImagePointsArray<double> > testPts;
testPts->InitializeArray(img.GetPointer());
testPts->SetName("pts");

vtkNew<vtkPoints> points;
points->SetData(testPts.GetPointer());

vtkNew<vtkCellSet> cellSet;
cellSet->SetPoints(points.GetPointer());
cellSet->GetPointData()->PassData(wavelet->GetPointData());
cellSet->GetImplementation()->SetDataSet(wavelet);

vtkNew<vtkIdTypeArray> ids;
vtkIdType ncells = wavelet->GetNumberOfCells();
ids->SetNumberOfTuples(ncells/2);
for (vtkIdType i=0; i<ncells/2; i++)
  {
  ids->SetValue(i, i);
  }

cellSet->GetImplementation()->SetCellIds(ids.GetPointer());

vtkNew<vtkContourFilter> contour;
contour->SetValue(0, 150);
contour->SetInputData(cellSet.GetPointer());
contour->Update();

Here, we restricted the contour filter to the first half of the cells of the original image data by simply creating a vtkCellSet. Note that this code would have been almost identical whether we were dealing with an image data, structured grid or an unstructured grid. All vtkCellSet needs is a vtkDataSet and a set of cell ids.

This is it! When I first worked on this example, I could not believe that it is this simple. Granted, I left some of the more complex methods such as GetPointCells() unimplemented because I didn't need them. Still, with some additional data structures, it would be fairly easy to implement them. We will discuss a new threshold filter that uses the vtkCellSet class in a future blog.