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

[dgobbi (David Gobbi)]

@zhang-qiang-github

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]?

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.

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

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

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?

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
]