We are trying to handle CT reconstruction with variable z spacing. In order to work with our software, we must obtain a vtkImageData with fixed spacing in all direction.
What we are doing for now is the following :
→ create a vtkRectilinearGrid from the dicom slices
→ use vtkResampleToImage to interpolate the grid and create a vtkImageData.
The problem is that vtkResampleToImage is really slow, in order to interpolate a grid of 512x512x556 into and image of 512x512x556 it takes up to 40 seconds on an intel core i7-9700K (Windows 64, vtk 9.0.1)
Is there any other filter that can interpolate a rectilinear grid into an image data ? Or is there some option that can be set to accelerate the interpolation (the filter seems to be using vtkProbeFilter under the hood) ?
And a final question, we also use ITK and OpenCV, does anyone knows if there is a faster filter to achieve this in those libraries ?
@simonesneault It might help to understand what is your use case for the data after resampling. If you’re just rendering it (as a volume, say), vtk can render rectilinear grids without the need to resample. There are other filters that can do slicing, etc directly on the rectilinear grid.
ITK has itk::BSplineScatteredDataPointSetToImageFilter, which is a very generic filter for resampling scattered data to create an image. I can’t say anything about the speed, however.
Hello Sankhesh,
Thanks for your answer
Well the use case is simple, some of our clients are now sending us some CTs with variable inter slice (it seems to become not so unusual), and we need a vtkImageData because we do rendering and slicing (MPR) but also segmentation and processing using ITK. Our 12+ years softwares are entirely based on the structure of the vtkImageData/itkImage and changing to a vtkRectilinearGrid seems quite a challenge, if not impossible … I don’t think that ITK filters can work with variable spacing, but maybe I’m wrong ?
@pieper implemented an amazingly simple and very fast solution for this in 3D Slicer: create a grid transform to transfer the 4 corners of each slice into the correct physical location. You don’t even have to apply this to the entire volume (unless you want to use volume rendering), but you can use this transform to get slices as needed, using vtkImageReslice filter. You can get the correct slice (or volume) if slices are unevenly spaced, slices are missing, volume is sheared, and even if slices are not all parallel (as long as slices not intersect, as intersection would prevent computation of the inverse transform).
Hello @lassoan
I’ve just checked with 3D slicer and it can load properly CT’s with variable spacing. There don’t seem to be any interpolation at all, but i’ve read the comment on the relevant PR here and understand the motivations/argumentations for that. I might end up with that kind of simple solution, ie replicate the nearest slice without interpolation …
Thanks for the advice
@lassoan’s suggestion of using @pieper’s solution for transforming just the grid corners is an excellent idea.
For interpolation, an alternative solution I can think of is using the volume mapper to render the rectilinear grid with blend mode set to slice_blend. You can set enable the mapper’s render to texture mode which would give you vtkImageData for the slice.
The grid transform and image reslice filter performs linear interpolation by default. The volume raycast mappers use linear texture interpolation, too, so I don’t think it would make a difference.
If needed, you can easily switch to cubic or spline interpolation in the reslice filter. The image quality could be further improved in some cases by using bspline transform instead of grid transform, but that might introduce complexity and unexpected behaviour, so I’m not sure if it is worth going into that.
The default interpolation is nearest-neighbor, not linear.
The vtkImageBSplineInterpolator must only be used on images that have been pre-filtered with the vtkImageBSplineCoefficients filter. Otherwise the result will be blurry.
Likewise, vtkBSplineTransform cannot be used directly with a displacement grid. The displacement grid must be filtered first (with vtkImageBSplineCoefficients). But for fixing missing slices, I’m pretty sure that a linearly-interpolated grid (with vtkGridTransform) is what you want.
!