I have some data (electron density) in a structured grid. I can plot the isosurface for different isovalues, and it looks nice. I’d like to integrate the density inside (or outside) the isosurface. Does VTK offer some tool for that? And would it be more accurate than a plain “rectangle rule” integration (sum of all values larger/smaller than the isovalue, multiplied by the volume element of the grid)? Any caveats if, for example, the isosurface is not close (because it’s clipped at the grid limits)?
Hi @Jellby,
Yes! there is a filter for that - vtkIntegrateAttributes
It looks like you’re interested in the volume integral of the electron density? One possible approach is you can clip your structured grid with a “scalar implicit function = isovalue” and pass the clipped volume to vtkIntegrateAttributes.
I think it should be fine because the cells at grid limits would still be whole (like a voxel). As long as you do not end up with lower dimensional cells after clipping, it should be okay.
Thanks! There seems to be some lack of examples for the usage of vtkIntegrateAttributes, so this is what I used:
clip = vtk.vtkClipDataSet()
clip.SetInputConnection(xyz.GetOutputPort())
integ = vtk.vtkIntegrateAttributes()
integ.SetInputConnection(clip.GetOutputPort())
# later, when I want to do the integration, and "clip" could have gone out of scope, but not "integ":
integ.GetInputConnection(0, 0).GetProducer().SetValue(isovalue)
integ.Update()
print(integ.GetOutput().GetPointData().GetArray(0).GetTuple(0)[0])
However, the numbers I get are very close to the naive integration:
vol = np.prod(xyz.GetInput().GetSpacing())*abs(xyz.GetTransform().GetMatrix().Determinant())
data = numpy_support.vtk_to_numpy(xyz.GetInput().GetPointData().GetScalars())
print(np.sum(data[data >= isovalue])*vol)
so I’m not sure they’re “better” (or faster).
Hmm. If you’re looking for accuracy, piecewise quadratic integration would be the best. I do not think Simpson’s rule is implemented anywhere in VTK. It would be a great addition to that filter.