Converting Pointer-Based Algorithms to Zero Copy Arrays

In my last blog, I demonstrated how VTK’s new “zero copy” infrastructure can be used to develop data arrays that use memory layouts different than VTK’s default. Many of VTK’s algorithm will work out of box with these arrays – the ones that depend on virtual methods such as SetTuple()/GetTuple() and SetValue()/GetValue(). However, certain algorithms need to be modified slightly. We’ll cover how this can be achieved here.

If you pass a subclass of vtkMappedDataArray<> to an algorithm that accesses the raw data pointer of an array via the GetVoidPointer() method or its ilk, you will see the following error at runtime:

Warning: In /Users/berk/Work/VTK/git/Common/Core/vtkMappedDataArray.txx, line 47
19vtkImagePointsArrayIdE (0x7fa939c0ecd0): GetVoidPointer called.
This is very expensive for vtkMappedDataArray subclasses, since the
scalar array must be generated for each call. Consider using a
vtkTypedDataArrayIterator instead.

If you see this error, it may be a good idea to refactor the offending algorithm to avoid this deep copy (which is slow and wastes memory). As an example, let’s consider a piece of code that uses the vtkImagePointsArray from the last blog.

vtkNew<vtkRTAnalyticSource> source;
source->SetWholeExtent(-100, 100, -100, 100, -100, 100);
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<vtkStructuredGrid> sgrid;
sgrid->SetDimensions(wavelet->GetDimensions());
sgrid->SetPoints(points.GetPointer());
sgrid->GetPointData()->SetScalars(wavelet->GetPointData()->GetScalars());

vtkNew<vtkGridSynchronizedTemplates3D> contour2;
contour2->SetInputData(sgrid.GetPointer());
contour2->SetValue(0, 200);
contour2->Update();

The idea behind this is to treat an image data as a vtkStructuredGrid without creating an explicit point array. Not very useful on its own but it can be extended fairly easily to represent dataset not currently efficiently handled by VTK such as a grid curvilinear in x-y but regular in z, a regular spherical grid etc. When you run this example, you will notice the warning I mentioned above. The result will be correct, however. If you track down the origin of the warning, you end up in vtkGridSynchronizedTemplates3D() :

template <class T, class PointsType>
void ContourGrid(vtkGridSynchronizedTemplates3D *self, ...)
{
  int *inExt = input->GetExtent();
  int xdim = exExt[1] - exExt[0] + 1;
  int ydim = exExt[3] - exExt[2] + 1;
  double n0[3], n1[3];  // used in gradient macro
  double *values = self->GetValues();
  int numContours = self->GetNumberOfContours();
  PointsType *inPtPtrX, *inPtPtrY, *inPtPtrZ;
  PointsType *p0, *p1, *p2, *p3;
  T *inPtrX, *inPtrY, *inPtrZ;
  T *s0, *s1, *s2, *s3;
  int XMin, XMax, YMin, YMax, ZMin, ZMax;
  int incY, incZ;
  PointsType* points =
    static_cast<PointsType*>(
    input->GetPoints()->GetData()->GetVoidPointer(0));

The problem is in the GetVoidPointer(0) call that asks for a raw pointer of type PointsType (double in this case). The solution to this issue is described in detail in David Lonie’s excellent Wiki page. Using the method described by David, I converted vtkGridSynchronizedTemplates3D. The full code can be found here. The change is minimal. I first made a change to the following function:

template <class T>
void ContourGrid(vtkGridSynchronizedTemplates3D *self,
                 int *exExt, T *scalars, vtkStructuredGrid *input,
                 vtkPolyData *output, vtkDataArray *inScalars,
                 bool outputTriangles)
{
  switch(input->GetPoints()->GetData()->GetDataType())
    {
    vtkTemplateMacro(
      ContourGrid(self, exExt, scalars, input, output,
        static_cast<VTK_TT *>(0), inScalars, outputTriangles));
    }
}

For those not familiar with VTK’s template macros, vtkTemplateMacro expands to a switch statement that looks like the following.

#define vtkTemplateMacroCase(typeN, type, call)     \
  case typeN: { typedef type VTK_TT; call; }; break
#define vtkTemplateMacro(call)                                              \
  vtkTemplateMacroCase(VTK_DOUBLE, double, call);                           \
  vtkTemplateMacroCase(VTK_FLOAT, float, call);                             \
  vtkTemplateMacroCase_ll(VTK_LONG_LONG, long long, call)                   \
  vtkTemplateMacroCase_ll(VTK_UNSIGNED_LONG_LONG, unsigned long long, call) \
...

This allows calling a templated function such as ContourGrid with the right template type. See above for the original ContourGrid signature with 2 template arguments. I changed the single template argument ContourGrid to the following.

template <class T>
void ContourGrid(vtkGridSynchronizedTemplates3D2 *self,
                 int *exExt, T *scalars, vtkStructuredGrid *input,
                 vtkPolyData *output, vtkDataArray *inScalars,
                 bool outputTriangles)
{
  vtkDataArray* pts = input->GetPoints()->GetData();
  switch(pts->GetDataType())
    {
    vtkDataArrayIteratorMacro(pts,
      ContourGrid(self, exExt, scalars,
        input, output, static_cast<vtkDAValueType*>(0),
        inScalars, outputTriangles, vtkDABegin));
    }
}

Note the usage of vtkDataArrayIteratorMacro macro. This macro was introduced by David along with vtkDABegin and vtkDAValueType. See his Wiki page for details. I then changed the other ContourGrid method as follows.

template <class T, class PointsType, class InputIterator>
void ContourGrid(vtkGridSynchronizedTemplates3D2 *self,
                 int *exExt, T *scalars,
                 vtkStructuredGrid *input, vtkPolyData *output,
                 PointsType*, vtkDataArray *inScalars, bool outputTriangles,
                 InputIterator points)
{
...
}

Note how points is now of type InputIterator which is a template argument magically defined within the vtkDataArrayIteratorMacro. This class leverages the vtkMappedDataArray<> API to provide pointer style semantics. This is it! After these changes, the pipeline shown above can now be run without any deep copies.

In future blogs, I will continue to expand on zero copy structures as well as dig into the cost of using these abstractions in terms of performance.