Select columns of structured data

ParaView Support
1

I’m working with image data in a programmable filter, and I would like to work with entire columns or sheets of celldata (e.g. data[inda, indb, :], data[inda, :, :]) Is there a way for me to do this?

2

Hi Nathan,

Could you please clarify your question or provide some data? If you want to use all the values of one cell property, you can get them into numpy array. If you want to use all the values of all the cell properties, you can combine them in 2d array. Maybe, I do not understand the question.

Best Regards,
Pavel

3

What I want to do is manipulate a single CellData field as if it were a multidimensional, dense, array. It sounds like I can’t do that, though. It wouldn’t be possible with any data type other than a uniform rectilinear grid, anyway.

I expect I can use indexing to work around this, though? Something like:

nx, ny, nz = field.GetDimensions()
for indx in range(nx):
  for indy in range(ny):
    new_field[indx, indy] = sum(field[indx*nx*ny + indy*ny:indx*nx*ny + indy*ny + nz])

I’m not thrilled about the idea of working out the details of that indexing strategy, though. Maybe it’s as simple as:
3d_field = ndarray((nx,ny,nz), buffer = field, order='C')
These approaches are making a lot of assumptions about how vtkImageData is laid out in memory, though.

4

Thanks for the clarification. If I undestand you correctly, you can use vtk_to_numpy python module from vtk.util.numpy_support to avoid this problem of thinking about “how vtkImageData is laid out in memory”.

5

E.g. if img is a an object of vtkDataImage, you could use something like this:

from vtk.util.numpy_support import vtk_to_numpy

nx, ny, _ = img.GetDimensions()
field_vtk = img.GetCellData().GetScalars()
field_numpy = vtk_to_numpy(field_vtk)
3d_field_numpy = field_numpy.reshape(nx, ny, -1)

I’m not sure there are no errors in a code above (I don’t have an example dataset), but hope it helps.

6

So this looks like it’s almost working (as a programmable filter), but the end result is about half what it should be. I’m trying to integrate a 3D density field with respect to ‘Z’.

I can’t for the life of me figure out why I’d need a factor of 2, so I’m wondering if I’ve messed up the data structure interface somehow.

from vtk.util.numpy_support import vtk_to_numpy, numpy_to_vtk
nx, ny, nz = inputs[0].GetDimensions()
dx, dy, dz = inputs[0].GetSpacing()
integrand = vtk_to_numpy(inputs[0].PointData['density']).reshape(nz, ny, nx)
out_numpy = numpy.add.reduce(integrand, axis=0, keepdims=True)*dz
out_vtk = numpy_to_vtk(out_numpy.flatten(), deep=1)
out_vtk.SetName('density_a')
out_vtk.SetNumberOfComponents(1)
out_vtk.SetNumberOfTuples(nx*ny*1)
out = self.GetOutput()
out.SetDimensions(nx, ny, 1)
out.SetOrigin(inputs[0].GetOrigin())
out.SetSpacing(inputs[0].GetSpacing()) 
out.GetPointData().AddArray(out_vtk)
1 Like
7

Could you please clarify:

Could you share some data?

Edited: code looks pretty much clear. Could you please compare the sum of all initial point values (vtk_to_numpy(inputs[0].PointData[‘density’]) multiplied by dz with the sum of all elements of out_vtk array in the end? Are they the same or is factor of 2 needed?

Edited 2: code from your previous message works fine for my data. Integration is ok. No need for factor of 2.

1 Like