[vtkusers] vtkImagePlaneWidget interpolation problem

[David Gobbi]

Hi Oliver,

The extra squares are due to an extra level of interpolation, you can get
rid of them by using this:

    plane_widget.SetResliceInterpolateToNearestNeighbour()

Also, in your code you have assumed that the origin of the image data
volume is the outer corner of the first voxel.  This is not correct.  In
VTK, the origin is always the center of the first voxel.  As a result, you
have placed the plane at the wrong position.  Use this instead:

    plane_widget.SetOrigin(-0.5, -0.5, 0.0)
    plane_widget.SetPoint1(+7.5, -0.5, 0.0)
    plane_widget.SetPoint2(-0.5, +7.5, 0.0)

    create_sphere(-0.5, -0.5, 0, (1, 0, 0))
    create_sphere(+2.5, -0.5, 0, (0, 1, 0))
    create_sphere(-0.5, +2.5, 0, (0, 0, 1))
    create_sphere(+2.5, +2.5, 0, (1, 1, 0))

Note that DICOM, NIfTI, and other medical image formats also measure
positions from the center of the first voxel.


Now, here is some info about the vtkImagePlaneWidget.  It's an old class
and is a little clunky.  It uses two steps:
1) first, vtkImageReslice extracts a slice (often padded or cropped) from
the volume
2) next, the output of vtkImageReslice is used as a texture for a
rectangular polygon

Both the reslicing and the texture mapping perform interpolation (the
interpolation method for either can be controlled).

The main problem with vtkImagePlaneWidget is that the texels don't always
match the image voxels.  In fact, they only match under very specific
circumstances:

1) the corner of the plane must be at the outer corner of the first voxel
(i.e. a half-voxel offset in X and Y from the volume origin, or a multiple
of the voxel spacing from the aforementioned position)

2) the size of the plane must be a power of two of the voxel size, because
this class was written in a time when GPUs had a power of two restriction
on texture sizes

With these restrictions in mind, if you have a volume made of slices that
are power-of-two (e.g. 256x256 or 512x512), then as long as the plane is
placed as specified above, it is possible to exactly render the image
slices with vtkImagePlaneWidget.


I've written a mapper called vtkImageResliceMapper that overcomes the voxel
vs. texel issue by putting the texels at exactly the same positions as the
voxels (it also has another mode where the texels are placed exactly at
screen pixel locations), but as far as I know no-one has used this class to
make a VTK widget.

 - David




On Thu, Sep 28, 2017 at 1:31 PM, Oliver Natt <onatt at gmx.de> wrote:

> Dear vtkusers,
>
> during my work on a small programm for the analysis of computed tomography
> data I found an unexpected behaviour of the vtkImagePlaneWidget. I tried to
> boil down the problem to a small python script, which uses a very small
> image dimension of only 3 x 3 x 3, where the effect becomes prominent (see
> python program below).
>
> The script makes a cut through a 3 x 3 x 3 volume using a
> vtkImagePlaneWidget. The four coloured spheres indicate the expected edges
> of the resulting image slice. The result is pretty much, what I expected
> (see attachement screenshot1.png). However, when I grab the
> vtkImagePlaneWidget using "Control + Middle Mouse Button" an move the plane
> around in the x-y plane, some unexpected effects occur (see
> screenshot2.png): Sometimes some kind of interpolation happens, which
> results in an image that looks rather like a 6 x 3 grid instead of 3 x 3
> and also the position of the image relative to the four coloured spheres
> changes. What I require is a cut through the volume data that preserves the
> correct voxel boundaries. Is that possible using the vtkImagePlaneWidget or
> is there an alternative class, that serves this purpose?
>
> Many thanks
> Oliver
>
> import vtkimport vtk.util.numpy_supportimport numpy as np
> # create a rendererrenderer = vtk.vtkRenderer()
> renderer.SetBackground(0.2, 0.2, 0.5)
> renderer.SetBackground2(0.4, 0.4, 1.0)
> renderer.SetGradientBackground(True)
> # create a render_windowrender_window = vtk.vtkRenderWindow()
> render_window.AddRenderer(renderer)
> # create a renderwindowinteractoriren = vtk.vtkRenderWindowInteractor()
> iren.SetRenderWindow(render_window)
> iren.Initialize()
> # Create a 3 x 3 x 3 image data with a spacing of 1.image = vtk.vtkImageData()
> dat = np.linspace(1, 27, 27, dtype=np.int16) * 10dat = np.ascontiguousarray(dat)
> dat = vtk.util.numpy_support.numpy_to_vtk(dat, deep=1)
> image.SetDimensions(3, 3, 3)
> image.SetSpacing(1, 1, 1)
> image.GetPointData().SetScalars(dat)
> image.Modified()
> # Create a vtkImagePlaneWidget and activate itplane_widget = vtk.vtkImagePlaneWidget()
> plane_widget.TextureInterpolateOff()
> plane_widget.SetInputData(image)
> plane_widget.SetInteractor(render_window.GetInteractor())
> plane_widget.SetOrigin(0, 0, 0)
> plane_widget.SetPoint1(8, 0, 0)
> plane_widget.SetPoint2(0, 8, 0)
> plane_widget.UpdatePlacement()
> plane_widget.DisplayTextOn()
> plane_widget.On()
> plane_widget.SetWindowLevel(200, 50)
> rs = plane_widget.GetReslice()
>
> def create_sphere(x, y, z, col):
>     """Create a small sphere with color col at postion x, y, z"""    source = vtk.vtkSphereSource()
>     source.SetCenter(x, y, z)
>     source.SetRadius(0.2)
>     mapper = vtk.vtkPolyDataMapper()
>     mapper.SetInputConnection(source.GetOutputPort())
>     actor = vtk.vtkActor()
>     actor.SetMapper(mapper)
>     actor.GetProperty().SetColor(col)
>     renderer.AddActor(actor)
>
> # Position 4 small markers at the expected edges of the image slice.create_sphere(0, 0, 0, (1, 0, 0))
> create_sphere(3, 0, 0, (0, 1, 0))
> create_sphere(0, 3, 0, (0, 0, 1))
> create_sphere(3, 3, 0, (1, 1, 0))
> # Reset the camera and start the RenderWindowInteractorrenderer.ResetCamera()
> iren.Start()
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20170928/cbf69303/attachment.html>