How to convert the coordinate between pixel and physical (vtkImageReslice)?

Support
1

I reslice a 2D image from a 3D volume image by vtkImageReslice. Then, I convert the image (matrix， numpy) from vtkImageReslice.GetOutput(). My question is: for one pixel (x, y) in image, what’s the coordinate in 3D volume?

I have found a solution to convert the coordinate from vtkImageReslice to 3D volume: vtkImageReslice.GetResliceAxes().MultiplyPoint().

However, how to convert the coordinate from pixel to vtkImageReslice? The vtkImageReslice need two inputs: the axes and center point. However, I find that some different center point would have the same output. And my test code is:

from vtk.util.numpy_support import vtk_to_numpy, numpy_to_vtk 
import vtk 
import numpy as np 

def vtkToNumpy(data): 
    temp = vtk_to_numpy(data.GetPointData().GetScalars()) 
    dims = data.GetDimensions() 
    numpy_data = temp.reshape(dims[2], dims[1], dims[0]) 
    numpy_data = numpy_data.transpose(2,1,0) 
    return numpy_data 

def numpyToVTK(data): 
    flat_data_array = data.transpose(2,1,0).flatten() 
    vtk_data_array = numpy_to_vtk(flat_data_array) 
    vtk_data = numpy_to_vtk(num_array=vtk_data_array, deep=True, array_type=vtk.VTK_FLOAT) 
    img = vtk.vtkImageData() 
    img.GetPointData().SetScalars(vtk_data) 
    img.SetDimensions(data.shape) 
    return img 

img = np.zeros(shape=[512,512,120]) 
img[0:300,0:100,:] = 1 

vtkImg = numpyToVTK(img) 

reslice = vtk.vtkImageReslice() 
reslice.SetInputData(vtkImg) 
reslice.SetAutoCropOutput(True) 
reslice.SetOutputDimensionality(2) 
reslice.SetInterpolationModeToCubic() 
reslice.SetSlabNumberOfSlices(1) 
reslice.SetOutputSpacing(1.0,1.0,1.0) 

axialElement = [ 
    1, 0, 0, 256, 
    0, 1, 0, 1000, 
    0, 0, 1, 100, 
    0, 0, 0, 1 
] 
resliceAxes = vtk.vtkMatrix4x4() 
resliceAxes.DeepCopy(axialElement) 
reslice.SetResliceAxes(resliceAxes) 
reslice.Update() 

reslicedImg = reslice.GetOutput() 
reslicedNpImg = vtkToNumpy(reslicedImg) 
import matplotlib.pyplot as plt 
plt.figure() 
plt.imshow(reslicedNpImg[:,:,0]) 
plt.show()

In the above code, the center point is (256,1000,100). Actually, this center point is outside of the original image (512,512,120), but the output is fine. And the center point (256, x, 100) would output the same image no matter what the x is. By the way, how to get the pixel index (x, y) for the center point in the resliced image.

So, for a pixel (x, y) in reslicedNpImg, how to obtain the coordinate in 3D image?

Thank you very much in advance!

1 Like
2

I think the piece you are missing is that SetAutoCropOutput(True) causes vtkImageReslice to choose an Origin for the output slice:

print("spacing %g %g %g" % reslicedImg.GetSpacing())
print("origin %g %g %g" % reslicedImg.GetOrigin())
print("extent %i %i %i %i %i %i" % reslicedImg.GetExtent())

spacing 1 1 1
origin -256 -1000 0
extent 0 511 0 511 0 0

To transform a coordinate from ‘output pixel index’ to ‘input voxel index’ you must:

  1. multiply by the output spacing (which is 1,1,1) and add the output origin
  2. call vtkImageReslice.GetResliceAxes().MultiplyPoint()
  3. subtract the input origin (which is 0,0,0) and divide by the input spacing (which is 1,1,1)

For your situation, one can get rid of the identity operations to simplify this:

  1. add the output origin to the pixel index
  2. call vtkImageReslice.GetResliceAxes().MultiplyPoint()
3

Thank you very much for your kindly answer. After adding the output origin to the pixel index and calculation with MultiplyPoint(), I still find something wrong. Please see my code:

from vtk.util.numpy_support import vtk_to_numpy, numpy_to_vtk 
import vtk 
import numpy as np 

def vtkToNumpy(data): 
    temp = vtk_to_numpy(data.GetPointData().GetScalars()) 
    dims = data.GetDimensions() 
    numpy_data = temp.reshape(dims[2], dims[1], dims[0]) 
    numpy_data = numpy_data.transpose(2,1,0) 
    return numpy_data 

def numpyToVTK(data): 
    flat_data_array = data.transpose(2,1,0).flatten() 
    vtk_data_array = numpy_to_vtk(flat_data_array) 
    vtk_data = numpy_to_vtk(num_array=vtk_data_array, deep=True, array_type=vtk.VTK_FLOAT) 
    img = vtk.vtkImageData() 
    img.GetPointData().SetScalars(vtk_data) 
    img.SetDimensions(data.shape) 
    return img 

img = np.ones(shape=[512,512,120])

vtkImg = numpyToVTK(img) 

reslice = vtk.vtkImageReslice() 
reslice.SetInputData(vtkImg) 
reslice.SetAutoCropOutput(True)
reslice.SetOutputDimensionality(2) 
reslice.SetInterpolationModeToCubic() 
reslice.SetSlabNumberOfSlices(1) 
reslice.SetOutputSpacing(1.0,1.0,1.0) 

import math
x = [math.cos(math.pi/4),math.sin(math.pi/4),0]
y = [-math.sin(math.pi/4),math.cos(math.pi/4),0]
z = [0,0,1]
center =[100, 254.5, 40]
axialElement = [
    x[0], y[0], z[0], center[0],
    x[1], y[1], z[1], center[1],
    x[2], y[2], z[2], center[2],
    0,    0,    0,    1
]
resliceAxes = vtk.vtkMatrix4x4()
resliceAxes.DeepCopy(axialElement)
reslice.SetResliceAxes(resliceAxes)
reslice.Update()

reslicedImg = reslice.GetOutput()
print(reslicedImg.GetOrigin())
print(reslicedImg.GetDimensions())

origin = reslicedImg.GetOrigin()
pixel = [0,0]
coor = reslice.GetResliceAxes().MultiplyPoint([pixel[0]+origin[0], pixel[1]+origin[1], origin[2], 1])
print(coor)


reslicedNpImg = vtkToNumpy(reslicedImg)
import matplotlib.pyplot as plt
plt.figure()
plt.imshow(reslicedNpImg[:,:,0])
plt.show()

In the above code, when the center point is [100, 254.5, 40], for the origin of resliced image, the coordinate in 3D volume is: [255.5 -255.5 40]. I can understand this result. However, when the center point is [100, 255.5 40], the coordinate in 3D volume for the origin of resliced image is [255.5, -255.5000152]. The analytical 3D coordinate of origin should be [255.5 -255.5]. I don’t know why the result is different only when center point is [100 255.5 40]?

Moreover, if I have a world coordinate, how can I calculate the pixel location in the resliced image? My code is:

get the world coordinate

origin = reslice.GetOutput().GetOrigin()
pixel = [0,0]
resliceAxes = reslice.GetResliceAxes()
coor = resliceAxes.MultiplyPoint([pixel[0]+origin[0], pixel[1]+origin[1], origin[2], 1])
print(coor)

# get the pixel location 
invertAxes = vtk.vtkMatrix4x4()
vtk.vtkMatrix4x4.Invert(resliceAxes,invertAxes)
location = invertAxes.MultiplyPoint([coor[0], coor[1], coor[2], 1])
print(location)
print([location[0]-origin[0], location[1]-origin[1], location[2]-origin[2]])

Actually, the [location[0]-origin[0], location[1]-origin[1], location[2]-origin[2]] should be [0,0,0], but the result is: [-1.21e-5, -7.41e-7, 0]. Is that right?

In addition, you said that add the output origin to the pixel index. Is that mean: [x+origin[0], y+origin[1], origin[2]], the z of pixel should be 0?

4

You should expect small errors in some results, because floating-point arithmetic has finite precision. However, for double the error should be much smaller than what you saw. I think that you found a bug. The VTK python wrappers seem to convert the numbers to float when vtkMatrix4x4.MultiplyPoint() is called. Because of this, only single-precision is used instead of double-precision. I will file a bug report so that this will be fixed.

For now, please use vtkMatrix4x4.MultiplyDoublePoint() instead of vtkMatrix4x4.MultiplyPoint(). This will force the use of double-precision.

Use vtkMatrix4x4.MultiplyDoublePoint() and the error will be much smaller. But it will not always be zero, because floating-point math is not exact.

Yes. If reslice.SetOutputDimensionality(2) is used, then the z of the output pixels is always zero. If reslice.SetOutputDimensionality(3) is used, then z can be anything.

I recommend that you try vtkImageReslice without SetAutoCropOutput(True) and without SetOutputDimensionality(2). If you do this, you can set the OutputOrigin to zero, which will simplify the transformations. The following code shows the necessary changes to the ResliceAxes matrix elements:

reslice = vtk.vtkImageReslice() 
reslice.SetInputData(vtkImg) 
#reslice.SetAutoCropOutput(True)
#reslice.SetOutputDimensionality(2) 
reslice.SetInterpolationModeToCubic() 
reslice.SetSlabNumberOfSlices(1) 
reslice.SetOutputSpacing(1.0,1.0,1.0) 
reslice.SetOutputOrigin(0.0,0.0,0.0)
# choose your own output dimensions
reslice.SetOutputExtent(0,511,0,511,0,0)

import math
angle = math.pi/4
x = [math.cos(angle),math.sin(angle),0.0]
y = [-math.sin(angle),math.cos(angle),0.0]
z = [0.0,0.0,1.0]
# the center of the slice in input coordinates
center = [255.5, 255.5, 40.0]
# the center of the slice in output coordinates
ss = [255.5, 255.5, 0.0]
axialElement = [
    x[0], y[0], z[0], center[0] - ss[0]*x[0] - ss[1]*y[0] - ss[2]*z[0],
    x[1], y[1], z[1], center[1] - ss[0]*x[1] - ss[1]*y[1] - ss[2]*z[1],
    x[2], y[2], z[2], center[2] - ss[0]*x[2] - ss[1]*y[2] - ss[2]*z[2],
    0.0,    0.0,    0.0,    1.0
]
5

Thank you very much. It really help me a lot.