[vtkusers] problem with vtkTransformToGrid

[(Jessie) Ting Guo]

I am trying to use vtkTransformToGrid to generate a grid file, but my
program didn't work.  There are three inputs: landmarks in source image,
landmarks in target image and source image. This program should generate a
grid file. I want to create a thinkplatespline transform according to the
landmarks in the source image and the target image, then use
vtkTransfromToGrid to generate a grid file containing the displacement
vectors. I am not sure if these vtk classes I am using can do what I want.
Could you guys please give me some advices?

import sys, string, array, fileinput
import os
import paths
from vtkMINCWriter import *
from vtkMNIXFMReader import *
from vtkMINCReader import *


try:
    FileName0 = sys.argv[1]
    FileName1 = sys.argv[2]
    InputImageName = sys.argv[3]
    OutputImageName = sys.argv[4]
except:
    print
    print 'Usage: python',
sys.argv[0],'<File0.tag>','<File1.tag>','<InputImageFile.mnc>','<Output.mnc>
'
    print
    sys.exit()

reader = vtkImageReader()
reader.SetFileName(InputImageName)
reader.ReleaseDataFlagOff()
reader.SetDataByteOrderToLittleEndian()
reader.SetDataExtent(0,255, 0,255, 0,247)
reader.SetDataSpacing(-1.17188,-1.17188,1)
reader.SetDataOrigin(154.414,151.514,-125.1)
#reader.SetDataMask 0x7fff
reader.Update

list_pos0 = []
list_pos1 = []

tagFile0 = open(FileName0,'r')
length0 = len(tagFile0.readlines())
for line0 in fileinput.input(FileName0):
    line0 = string.strip(line0)
    split0 = string.split(line0,'"')
    if len(split0) > 1:
        coords0 = string.split(split0[0]," ")
        x0=string.atof(coords0[0])
        y0=string.atof(coords0[1])
        z0=string.atof(coords0[2])
        point_pos0 = (x0,y0,z0)
        list_pos0.append(point_pos0)

tagFile1 = open(FileName1,'r')
length1 = len(tagFile1.readlines())
for line1 in fileinput.input(FileName1):
    line1 = string.strip(line1)
    split1 = string.split(line1,'"')
    if len(split1) > 1:
        coords1 = string.split(split1[0]," ")
        x1=string.atof(coords1[0])
        y1=string.atof(coords1[1])
        z1=string.atof(coords1[2])
        point_pos1 = (x1,y1,z1)
        list_pos1.append(point_pos1)

p0 = vtkPoints()
p1 = vtkPoints()
i = 0
for i in range (43):
    p1.InsertNextPoint(list_pos1[i])
    p0.InsertNextPoint(list_pos0[i])

transform = vtkThinPlateSplineTransform()
transform.SetSourceLandmarks(p0)
transform.SetTargetLandmarks(p1)
transform.SetBasisToR()

gridThinPlate = vtkTransformToGrid()
gridThinPlate.SetInput(transform)
gridThinPlate.SetGridExtent(0,255, 0,255, 0,247)
gridThinPlate.SetGridSpacing(-1.17188,-1.17188,1)
gridThinPlate.SetGridOrigin(154.414,151.514,-125.1)
gridThinPlate.SetGridScalarTypeToUnsignedChar()

gridTransform = vtkGridTransform()
gridTransform.SetDisplacementGrid(gridThinPlate.GetOutput())
gridTransform.SetDisplacementShift(gridThinPlate.GetDisplacementShift())
gridTransform.SetDisplacementScale(gridThinPlate.GetDisplacementScale())

reslice = vtkImageReslice()
reslice.SetInput(reader.GetOutput())
reslice.SetResliceTransform(gridTransform)
reslice.SetInterpolationModeToCubic()
reslice.SetOutputSpacing(1,1,1)
reslice.UpdateWholeExtent()

writer = vtkMINCWriter()
writer.SetFileName(OutputImageName)
writer.SetInput(reslice.GetOutput())
writer.SetOutputType('b')
writer.Write()