# How to convert the coordinate between pixel and physical (vtkImageReslice)?
**URL:** https://discourse.vtk.org/t/how-to-convert-the-coordinate-between-pixel-and-physical-vtkimagereslice/1391
**Category:** Support
**Created:** [July 21, 2019, 8:11am UTC](https://discourse.vtk.org/t/how-to-convert-the-coordinate-between-pixel-and-physical-vtkimagereslice/1391 "2019-07-21T08:11:40Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@zhang-qiang-github](https://discourse.vtk.org/u/zhang-qiang-github)
#### Post date: [July 21, 2019, 8:11am UTC](https://discourse.vtk.org/t/how-to-convert-the-coordinate-between-pixel-and-physical-vtkimagereslice/1391/1 "2019-07-21T08:11:40Z")
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:
```auto
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!
---
### Author:  [@dgobbi](https://discourse.vtk.org/u/dgobbi)
#### Post date: [July 22, 2019, 5:10pm UTC](https://discourse.vtk.org/t/how-to-convert-the-coordinate-between-pixel-and-physical-vtkimagereslice/1391/2 "2019-07-22T17:10:32Z")
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()
---
### Author:  [@zhang-qiang-github](https://discourse.vtk.org/u/zhang-qiang-github)
#### Post date: [July 23, 2019, 2:21am UTC](https://discourse.vtk.org/t/how-to-convert-the-coordinate-between-pixel-and-physical-vtkimagereslice/1391/3 "2019-07-23T02:21:52Z")
> [@zhang-qiang-github](#):
>
> 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()
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.2_1e-5, -7.4_1e-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?
---
### Author:  [@dgobbi](https://discourse.vtk.org/u/dgobbi)
#### Post date: [July 23, 2019, 7:43pm UTC](https://discourse.vtk.org/t/how-to-convert-the-coordinate-between-pixel-and-physical-vtkimagereslice/1391/4 "2019-07-23T19:43:32Z")
> [@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
]
```
---
### Author:  [@zhang-qiang-github](https://discourse.vtk.org/u/zhang-qiang-github)
#### Post date: [July 24, 2019, 2:43am UTC](https://discourse.vtk.org/t/how-to-convert-the-coordinate-between-pixel-and-physical-vtkimagereslice/1391/5 "2019-07-24T02:43:14Z")
> [@dgobbi](#):
>
> MultiplyDoublePoint
Thank you very much. It really help me a lot.