[Paraview] Finite-size voxels and Point vs. Cell Centered Image Data

[Eric Nodwell]

What is the standard way to represent image data in VTK/ParaView, where the image data is intended to represent finite sized voxels?  Here the 3D image data represent values at the centre of voxels (or presumably some sort of weighted value centred on the voxel centre).

This may seem like a trivial question, however I am unable to obtain the behaviour that seems most appropriate.  For example, suppose I want to display all voxels in a (3D) image that exceed a certain threshold value.  The obvious thing to do seems to be to save the data in a VTK .vti file and then use the Theshold filter in ParaView.  The Threshold filter will display Cells, or not display Cells, depending on the threshold.  But Cells of a vtkDataObject are hexagons connecting the points, not hexagons centred on the points.  (There is an option "All Scalars" for the Threshold filter which determines how the cell value is calculated from the neighbouring point values, but it has no effect on the location of the cell.)

One possibility may be that we are storing the image data incorrectly in the .vti file.  Right now we have are generating the data as vtkImageSource with Scalar Point Data.  Perhaps it would be more logical to generate and save this as vtkImageSource with Scalar Cell Data, and indeed the ParaView documentation seems to imply this (see http://paraview.org/OnlineHelpCurrent/XMLImageDataReader.html ).  However, vtkImageSource does not seem to be intended to be used in this way.  In the vtkImageData implementation, all the methods that access "scalar" data do not check for the existence of Scalars for Cell Data, but instead always access the Point Data Scalars, like this:

  vtkDataArray *scalars = this->GetPointData()->GetScalars();

Or to present this another way, consider the trivial example of a single voxel.  We would like to have this represented as a finite size box, but the data is inherently a single value.  Here is some example VTK python code that attempts to generate a single voxel:

  image = vtk.vtkImageData()
  image.SetScalarTypeToInt()
  image.SetSpacing (1, 1, 1)
  image.SetOrigin (0, 0, 0)
  image.SetDimensions (1, 1, 1)
  image.SetScalarComponentFromFloat (0,0,0,0,10)
  print image.GetBounds()

The output for the bounds is:

(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

This is not necessarily incorrect in the right context.  One example of such a context might be the case where the image data represent discrete samples taken at distinct points in space.  However it does seem inappropriate for data meant to represent voxels as small volume elements (e.g. most medical images), where we would like a single voxel to have finite real-space bounds and be displayed as such.  I have again the feeling that vtkImageData is not the appropriate object to use to represent voxel data, but I don't know what the correct alternative would be (something with scalar cell data presumably).

As an aside, if you save the above image object with vtkXMLImageDataWriter, the file can't even be opened in ParaView 3.18.0, as the following error occurs:

ERROR: In /Users/eric/source/ParaView-3.8.0/Servers/Filters/vtkTexturePainter.cxx, line 286
vtkTexturePainter (0x11de86350): Incorrect dimensionality.

However, a vtkImageData object with 2x2x2 dimensions can be opened and will display (using for example the Threshold filter) as a 1x1x1 box in ParaView.

Any insight would be appreciated.

Thanks,
Eric