View large images with vtkImageResliceMapper

ViewImage2D.py

#! /usr/bin/env python

import sys
import os
import mmap

import vtk

usage = """
Usage: "ViewImage2D [options] image
Options:
 -n   Nearest-neighbor interpolation
 -c   Cubic interpolation
 -s   Sinc interpolation (high-quality)
 -a   Antialias with band-limited sinc (slow)
 -r   Red channel only
 -g   Green channel only
 -b   Blue channel only
Image types:
  .png, .jpg, .tif, .lsm
Note: some tiff features might be unsupported.
"""

# count the number of filenames on the command line
filenames = []
for arg in sys.argv[1:]:
    if arg[0] != '-':
        filenames.append(arg)

if len(filenames) == 0:
    sys.stdout.write(usage)
    sys.exit(0)

name = filenames[0]
base,ext = os.path.splitext(name)

if ext in ('.jpg', '.jpeg'):
    reader = vtk.vtkJPEGReader()
elif ext in ('.tif', '.tiff', ".lsm"):
    reader = vtk.vtkTIFFReader()
elif ext in ('.png',):
    reader = vtk.vtkPNGReader()
elif ext in ('.bin',):
    # example for reading raw file with known header size and dims
    reader = vtk.vtkImageReader2()
    reader.SetHeaderSize(8394)
    reader.SetFileDimensionality(2)
    reader.SetDataScalarTypeToUnsignedShort()
    reader.SetDataByteOrderToLittleEndian()
    reader.SetDataExtent(0, 107, 0, 251, 0, 0)
    reader.SetDataSpacing(2.0, 1.0, 1.0)
else:
    sys.stderr.write("cannot read file type %s\n" % (ext,))
    sys.exit(1)

if len(filenames) == 1:
    reader.SetFileName(filenames[0])
    reader.Update()
else:
    filename_array = vtk.vtkStringArray()
    filename_array.SetNumberOfValues(len(filenames))
    for i,f in enumerate(filenames):
        filename_array.SetValue(i, f)
    reader.SetFileNames(filename_array)
    reader.Update()

"""
# example of memory-mapping a file into VTK
f = open("world.topo.bathy.200407.3x86400x43200.bin", "rb")
m = mmap.mmap(f.fileno(),0,mmap.MAP_PRIVATE,mmap.PROT_READ)
reader = vtk.vtkImageImport()
reader.SetImportVoidPointer(m)
reader.SetDataScalarTypeToUnsignedChar()
reader.SetNumberOfScalarComponents(3)
reader.SetDataExtent(0,86400-1,0,43200-1,0,0)
reader.SetWholeExtent(0,86400-1,0,43200-1,0,0)
reader.SetDataSpacing(1,1,1)
"""

# compute histogram (also forces entire image into memory)
hist = vtk.vtkImageHistogramStatistics()
hist.SetInputConnection(reader.GetOutputPort())
hist.Update()

# print some information about the image
size = reader.GetOutput().GetDimensions()
spacing = reader.GetOutput().GetSpacing()
if size[2] == 1:
    print("Size: %d %d" % size[0:2])
    print("Spacing: %f %f" % spacing[0:2])
else:
    print("Size: %d %d" % size)
    print("Spacing: %f %f" % spacing)

centerer = vtk.vtkImageChangeInformation()
centerer.SetInputConnection(reader.GetOutputPort())
centerer.CenterImageOn()

# for playing with the number of threads used
#vtk.vtkMultiThreader.SetGlobalMaximumNumberOfThreads(1)

# create the RenderWindow and Renderer
ren1 = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren1)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)

# set up the interpolator
interp = vtk.vtkImageInterpolator()
interp.SetInterpolationModeToLinear()

if '-s' in sys.argv:
    interp = vtk.vtkImageSincInterpolator()
    interp.SetWindowFunctionToKaiser()
    interp.SlidingWindowOn()
if '-a' in sys.argv:
    interp = vtk.vtkImageSincInterpolator()
    interp.SetWindowFunctionToKaiser()
    interp.AntialiasingOn()
    interp.SlidingWindowOn()
if '-r' in sys.argv:
    interp.SetComponentOffset(0)
    interp.SetComponentCount(1)
if '-g' in sys.argv:
    interp.SetComponentOffset(1)
    interp.SetComponentCount(1)
if '-b' in sys.argv:
    interp.SetComponentOffset(2)
    interp.SetComponentCount(1)

# create the mapper, attach the interpolator
im = vtk.vtkImageResliceMapper()
im.SetInterpolator(interp)
im.SliceFacesCameraOn()
im.SliceAtFocalPointOn()
im.JumpToNearestSliceOn()
#im.SeparateWindowLevelOperationOff()
im.ResampleToScreenPixelsOn()
im.SetInputConnection(centerer.GetOutputPort())
im.BorderOn()

# create image property manager
# (interpolation set here will not override mapper's interpolator)
ip = vtk.vtkImageProperty()
ip.SetInterpolationTypeToLinear()
if '-n' in sys.argv:
    ip.SetInterpolationTypeToNearest()
if '-c' in sys.argv:
    ip.SetInterpolationTypeToCubic()

# automatically set the window/level from the histogram
r = hist.GetAutoRange()
if r[1] - r[0] > 255:
    ip.SetColorWindow(r[1]-r[0])
    ip.SetColorLevel(0.5*(r[0]+r[1]))

# the vtkImageSlice is like vtkActor, but for images
ia = vtk.vtkImageSlice()
ia.SetMapper(im)
ia.SetProperty(ip)

# window size, increase as desired
renWin.SetSize(1200, 800)

# provide basic interaction with image
iren = vtk.vtkRenderWindowInteractor()
style = vtk.vtkInteractorStyleImage()
style.SetInteractionModeToImage2D()
if size[2] > 1:
    style.SetInteractionModeToImageSlicing()
iren.SetInteractorStyle(style)
renWin.SetInteractor(iren)
#renWin.FullScreenOn()

# use AddObserver to add our own event handler
def event_handler(style, event):
    if event == "LeftButtonPressEvent":
        style.StartPan()
    elif event == "MiddleButtonPressEvent":
        if size[2] > 1:
            style.StartSlice()
        else:
            pass
    elif event == "RightButtonPressEvent":
        pass

style.AddObserver("LeftButtonPressEvent", event_handler)
style.AddObserver("MiddleButtonPressEvent", event_handler)
#style.AddObserver("RightButtonPressEvent", event_handler)

# can use slice thickness for camera clipping range
if size[2] > 1:
    thickness = spacing[2]*(size[2]-1)
else:
    thickness = 1.0

# render the image
renWin.Render()
rsize = renWin.GetSize()
cam1 = ren1.GetActiveCamera()
cam1.ParallelProjectionOn()
cam1.SetParallelScale(0.5*spacing[1]*rsize[1])
cam1.SetFocalPoint(0.0, 0.0, 0.0)
cam1.SetPosition(0.0, 0.0, thickness)
cam1.SetClippingRange(0.5*thickness, 1.5*thickness)

ren1.AddViewProp(ia)
ren1.SetBackground(0.1,0.2,0.4)

renWin.Render()
iren.Start()

print("done!")