[vtkusers] Resizing dicom files to same volume

[Wright James (NHS Greater Glasgow & Clyde)]

VTK'ers

I have multiple DICOM datasets with different spacing and dimensions. I am using vtkImagePlaneWidget() to view them

 I want to resample the datasets so they occupy the same volume in space, to then use e.g. vtkImageBlend to merge the images by opacity. I have tried vtkImageResample and vtkImageReslice (in the documentation for this class it even mentions the method SetInformationInput() as the way to go for matching voxel sampling. However whatever I try still results in the final output having image slices of different scales on the screen. Could anyone guide me as to what pipeline to use with  to get the datasets aligned?

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------


import vtk


def orthoPlane(imageData, orientation, dimension, volPicker, renderer,
               interactor, lut, opacity, displayText = 1, enabled = 1):

    '''
    a function to take care of rendering an image plane on the scene - we use
    the vtk provided vtkImagePlaneWidget to handle this
    '''

    ipwx = vtk.vtkImagePlaneWidget()
    ipwx.SetInput(imageData)
    ipwx.SetPicker(volPicker)
    ipwx.GetTexturePlaneProperty().SetOpacity(opacity)

    ipwx.GetColorMap().SetLookupTable(lut)
    ipwx.SetCurrentRenderer(renderer)
    ipwx.SetInteractor(interactor)
    ipwx.PlaceWidget()
    if orientation == "sagittal":
        ipwx.SetPlaneOrientationToXAxes()
    elif orientation == "coronal":
        ipwx.SetPlaneOrientationToYAxes()
    elif orientation == "axial":
        ipwx.SetPlaneOrientationToZAxes()
    ipwx.SetSliceIndex(dimension/2)
    ipwx.SetDisplayText(displayText)
    ipwx.SetEnabled(enabled)

    return ipwx


def MoveCallback(obj,event):

    '''
    function that forces the two modalities to move together from one user
    interaction
    '''

    coronalSPECT.SetSlicePosition(coronalCT.GetSlicePosition())
    coronalSPECT.UpdatePlacement()



''' the vtk pipeline from here '''
if __name__ == "__main__" :


    # read the data
    folder1 = "/home/jamie/VTK/vtkPython/vtk_python/brain_data/CT"
    reader1 = vtk.vtkDICOMImageReader()
    reader1.SetDirectoryName(folder1)
    reader1.Update()
    W1,H1,D1 = reader1.GetOutput().GetDimensions()

    a1,b1 = reader1.GetOutput().GetScalarRange()
    #print "Range of image (Houndsfield units): %d--%d" %(a1,b1)

    folder2 = "/home/jamie/VTK/vtkPython/vtk_python/brain_data/SPECT"
    reader2 = vtk.vtkDICOMImageReader()
    reader2.SetDirectoryName(folder2)
    #reader2.SetDataSpacing(0.59,0.59,2.5)
    #reader2.SetDataDimensions(512,512,61)
    reader2.Update()

    #### ---- resample SPECT data to CT dimensions here ---- ####

    #reader2.GetOutput().SetSpacing((0.59,0.59,2.5))
    #reader2.GetOutput().SetDimensions((512,512,61))
    #reader2.Update()

#    # resize the data
#    resize = vtk.vtkImageResize()
#    resize.SetInputConnection(reader2.GetOuputPort())
#    resize.SetResizeMethodToMagnificationFactors()
#    resize.SetMagnificationFactors(0.5,0.5,1)
#    resize.InterpolateOn()
#    resize.Update()
#
#    # reslice the PET data to line up with the CT volume
#    reslice = vtk.vtkImageReslice()
#    #reslice.SetInformationInput(reader1.GetOutput())
#    reslice.SetInputConnection(reader2.GetOutputPort())
#    #reslice.SetOutputDimensionality(2)
#    reslice.SetInterpolationModeToLinear()
    # resize the data
#    magnify = vtk.vtkImageResample()
#    magnify.SetDimensionality(3)
#    magnify.SetInput(reader2.GetOutput())
#    #magnify.SetAxisOutputSpacing(0,0.59)
#    #magnify.SetAxisOutputSpacing(1,0.59)
#    magnify.SetAxisMagnificationFactor(1,2)
#    magnify.SetAxisMagnificationFactor(0,3)
#    magnify.SetAxisMagnificationFactor(2,2)
#    #magnify.SetAxisOutputSpacing(2,2.5)
#    #magnify.ReleaseDataFlagOff()
#    magnify.Update()

    W2,H2,D2 = reader2.GetOutput().GetDimensions()
    a2,b2 = reader2.GetOutput().GetScalarRange()

    # create the color transfer functions for the datasets
    # greyscale table for the CT
    lutCT = vtk.vtkColorTransferFunction()
    lutCT.AddRGBPoint(a1 , 0.0, 0.0, 0.0)
    lutCT.AddRGBPoint(b1 , 1.0, 1.0, 1.0)
    # basic rgb colortable/lut for the PET data
    lutPET = vtk.vtkColorTransferFunction()
    lutPET.AddRGBPoint(a2,         0.0, 0.0, 0.0)
    lutPET.AddRGBPoint(a2+(b2-a2)/4, 0.0, 0.5, 0.5)
    lutPET.AddRGBPoint(a2+(b2-a2)/2, 0.0, 1.0, 0.0)
    lutPET.AddRGBPoint(b2-(b2-a2)/4, 0.5, 0.5, 0.0)
    lutPET.AddRGBPoint(b2,         1.0, 0.2, 0.0)

    # renderer and render window
    ren = vtk.vtkRenderer()
    ren.SetBackground(.2, .2, .2)
    renWin = vtk.vtkRenderWindow()
    renWin.SetSize( 512, 512 )
    renWin.AddRenderer( ren )

    # render window interactor
    iren = vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow( renWin )

    # interactor style
    style = vtk.vtkInteractorStyleTrackballCamera()
    iren.SetInteractorStyle(style)

    # mouse picker to select info from the data - will be displayed on image
    picker=vtk.vtkCellPicker()
    picker.SetTolerance(0.005)

    # coronal actors - do the same for the other planes
    coronalCT = orthoPlane(reader1.GetOutput(), "coronal", H1, picker, ren,
                            iren,lutPET, 0.65)

    coronalSPECT = orthoPlane(reader2.GetOutput(), "coronal", H2, picker, ren,
                           iren,lutPET, 0.25)



    # create an outline of the dataset
    outline = vtk.vtkOutlineFilter()
    outline.SetInput( reader1.GetOutput() )
    outlineMapper = vtk.vtkPolyDataMapper()
    outlineMapper.SetInput( outline.GetOutput() )
    outlineActor = vtk.vtkActor()
    outlineActor.SetMapper( outlineMapper )

    # add the outline actor to the scene and render it
    ren.AddActor( outlineActor )
    renWin.Render()

    # position camera to create a nice starting viewing angle. Note: the order
    # of operations on the camera is important, i.e. a Roll,Elevation,Flip
    # gives a different orientation if we change the sequence of operations
    # i.e do the Roll second
    camera = ren.GetActiveCamera()
    camera.Roll(-30)
    camera.Elevation(110)
    camera.SetViewUp((0,0,-1))

    # add an observer to the event of one of the planes, using a priority 1 so
    # that it overrides the standard event. We then want to use the event to
    # update the one that is not being interacted with with the origin, point1
    # and point2 coordinates of the moving one.
    coronalCT.GetInteractor().AddObserver("MouseMoveEvent",MoveCallback,1)

    # initialize, and start the interactor
    iren.Initialize()
    iren.Start()

********************************************************************************************************************

This message may contain confidential information. If you are not the intended recipient please inform the
sender that you have received the message in error before deleting it.
Please do not disclose, copy or distribute information in this e-mail or take any action in reliance on its contents:
to do so is strictly prohibited and may be unlawful.

Thank you for your co-operation.

NHSmail is the secure email and directory service available for all NHS staff in England and Scotland
NHSmail is approved for exchanging patient data and other sensitive information with NHSmail and GSi recipients
NHSmail provides an email address for your career in the NHS and can be accessed anywhere

********************************************************************************************************************