Volume display - Interpolate turned off by default in newest stable?

Interpolated image shows you the faithfully reconstructed continuous signal from the discrete samples. This reconstruction is lossless if sample-rate criterion is satisfied (see more details here). The criterion is not always satisfied - that’s when you see staircase artifacts in the reconstructed image (in both slice views, volume rendering, and reconstructed surfaces):

Interpolation is still about the best thing you can do, and definitely provides much more faithful reconstruction of the original signal than no interpolation (= nearest neighbor interpolation). Also, you can always make the sample-rate criterion satisfied by applying a low-pass filter on the data. That’s why we always apply surface smoothing when surface is reconstructed from from labelmap volume.

You should be able to see small details better when interpolation is used. For example, try to draw the dark separating line in this image when it is displayed with/without interpolation:

You need to squint your eye and guess where to draw when voxels are shown without interpolation, while it is quite straightforward on the interpolated image.

I agree that it can be useful to show the voxel grid during segmentation so that you know how fine details you can paint and where. Note that this is the grid of the binary labelmap representation of the segmentation and not the grid of the background scalar volume.

We could quite easily add this grid display option to the segmentation displayable manager (above a certain zoom factor). You can simulate it with this Python code snippet:

import SampleData
volumeNode = SampleData.SampleDataLogic().downloadCTACardio()

import numpy as np

ijkToRasVtk = vtk.vtkMatrix4x4()
volumeNode.GetIJKToRASMatrix(ijkToRasVtk)
ijkToRas = slicer.util.arrayFromVTKMatrix(ijkToRasVtk)

volumeExtent = np.array(volumeNode.GetImageData().GetExtent(), dtype=float)
volumeExtent[0] -= 0.5
volumeExtent[2] -= 0.5
volumeExtent[4] -= 0.5
volumeExtent[1] += 0.5
volumeExtent[3] += 0.5
volumeExtent[5] += 0.5

dimensions = volumeNode.GetImageData().GetDimensions()

pts = vtk.vtkPoints()
quads = vtk.vtkCellArray()

for axisIndex in range(3):
    if axisIndex == 0:
        pointsIjk = np.array([
            [volumeExtent[0], volumeExtent[2], volumeExtent[4], 1.0],
            [volumeExtent[0], volumeExtent[3], volumeExtent[4], 1.0],
            [volumeExtent[0], volumeExtent[3], volumeExtent[5], 1.0],
            [volumeExtent[0], volumeExtent[2], volumeExtent[5], 1.0]]).T
        offsetIjk = np.array([[1, 0, 0, 0], [1, 0, 0, 0], [1, 0, 0, 0], [1, 0, 0, 0]], dtype=float).T
    elif axisIndex == 1:
        pointsIjk = np.array([
            [volumeExtent[0], volumeExtent[2], volumeExtent[4], 1.0],
            [volumeExtent[1], volumeExtent[2], volumeExtent[4], 1.0],
            [volumeExtent[1], volumeExtent[2], volumeExtent[5], 1.0],
            [volumeExtent[0], volumeExtent[2], volumeExtent[5], 1.0]]).T
        offsetIjk = np.array([[0, 1, 0, 0], [0, 1, 0, 0], [0, 1, 0, 0], [0, 1, 0, 0]], dtype=float).T
    else:
        pointsIjk = np.array([
            [volumeExtent[0], volumeExtent[2], volumeExtent[4], 1.0],
            [volumeExtent[1], volumeExtent[2], volumeExtent[4], 1.0],
            [volumeExtent[1], volumeExtent[3], volumeExtent[4], 1.0],
            [volumeExtent[0], volumeExtent[3], volumeExtent[4], 1.0]]).T
        offsetIjk = np.array([[0, 0, 1, 0], [0, 0, 1, 0], [0, 0, 1, 0], [0, 0, 1, 0]], dtype=float).T
    pointsRas = np.dot(ijkToRas, pointsIjk)
    offsetRas = np.dot(ijkToRas, offsetIjk)
    for planeIndex in range(dimensions[axisIndex]+1):
        quads.InsertNextCell(4)
        for pointIndex in range(4):
            quads.InsertCellPoint(pts.InsertNextPoint(pointsRas[0][pointIndex], pointsRas[1][pointIndex], pointsRas[2][pointIndex]))
        pointsRas += offsetRas

gridPoly = vtk.vtkPolyData()
gridPoly.SetPoints(pts)
gridPoly.SetPolys(quads)

gridModelNode = slicer.modules.models.logic().AddModel(gridPoly)
gridModelDisplayNode = gridModelNode.GetDisplayNode()
gridModelDisplayNode.SetVisibility2D(True)
gridModelDisplayNode.SetOpacity(0.3)
gridModelDisplayNode.SetVisibility3D(False)

# Running this line of code will show/hide the grid:
# gridModelDisplayNode.SetVisibility(not gridModelDisplayNode.GetVisibility())

See how the grid shows more relevant information than the non-interpolated master volume, especially when resolution of segmentation differs from the resolution of the master volume:

It is interesting to note that this way of showing the grid works well for arbitrarily reformatted slices:


It would be great if you could play around with the Python code snippet above and if you find grid display useful then we can discuss about adding it to segmentation displayable manager (similarly to how it worked in Slicer3).

1 Like