vtk: how to obtain the image pixel index from a world point

Question

If I pick a world point from a image, How can I convert the world coordinate to image index?

import vtk
import numpy as np
from vtk.util.numpy_support import numpy_to_vtk
def numpyToVTK(data, multi_component=False, type='float'):
    if type == 'float':
        data_type = vtk.VTK_FLOAT
    elif type == 'char':
        data_type = vtk.VTK_UNSIGNED_CHAR
    else:
        raise RuntimeError('unknown type')
    if multi_component == False:
        if len(data.shape) == 2:
            data = data[:, :, np.newaxis]
        flat_data_array = data.transpose(2,1,0).flatten()
        vtk_data = numpy_to_vtk(num_array=flat_data_array, deep=True, array_type=data_type)
        shape = data.shape
    else:
        assert len(data.shape) == 3, 'only test for 2D RGB'
        flat_data_array = data.transpose(1, 0, 2)
        flat_data_array = np.reshape(flat_data_array, newshape=[-1, data.shape[2]])
        vtk_data = numpy_to_vtk(num_array=flat_data_array, deep=True, array_type=data_type)
        shape = [data.shape[0], data.shape[1], 1]
    img = vtk.vtkImageData()
    img.GetPointData().SetScalars(vtk_data)
    img.SetDimensions(shape[0], shape[1], shape[2])
    return img
global sphereActor, textActor
sphereActor = None
textActor = None
def mouseMoveEvent(iren, event):
    x, y = iren.GetEventPosition()

    picker = vtk.vtkWorldPointPicker()
    picker.Pick(x, y, 0, render)
    worldPoint = picker.GetPickPosition()
    
    ##############################################
    ## convert world point to image index
    ##############################################

    sphere = vtk.vtkSphereSource()
    sphere.SetCenter(worldPoint[0], worldPoint[1], worldPoint[2])
    sphere.SetRadius(2)
    sphere.Update()
    sphereMapper = vtk.vtkPolyDataMapper()
    sphereMapper.SetInputData(sphere.GetOutput())
    global sphereActor, textActor
    if sphereActor != None:
        render.RemoveActor(sphereActor)
    sphereActor = vtk.vtkActor()
    sphereActor.SetMapper(sphereMapper)
    sphereActor.GetProperty().SetColor(255, 0, 0)
    render.AddActor(sphereActor)
    render.Render()

    if textActor != None:
        render.RemoveActor(textActor)
    textActor = vtk.vtkTextActor()
    textActor.SetInput('world coordinate: (%.2f, %.2f, %.2f)'%(worldPoint[0], worldPoint[1], worldPoint[2]))
    textActor.GetTextProperty().SetColor(1, 0, 0)
    textActor.GetTextProperty().SetFontSize(15)
    render.AddActor(textActor)


img = np.zeros(shape=[128, 128])
for i in range(128):
    for j in range(128):
        img[i, j] = i+j

vtkImg = numpyToVTK(img)
imgActor = vtk.vtkImageActor()
imgActor.SetInputData(vtkImg)

render = vtk.vtkRenderer()
render.AddActor(imgActor)
# render.Render()

renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(render)
renWin.Render()
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
iren.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera())
iren.Initialize()
iren.AddObserver('MouseMoveEvent', mouseMoveEvent)
iren.Start()

In the above code, if I don't rotate the image, the world point is (x, y, 0):

And it is agree with what I know. For the world point (x, y, z) and the image index (i, j, k), the conversion should be:

worldPoint (x,y,z) = i*spacingX*directionX + j*spacingY*directionY + k*spacingZ*directionZ + originPoint

In the above code, the image is converted from numpy, thus:

directionX = [1, 0, 0]
directionY = [0, 1, 0]
directionZ = [0, 0, 1]
originPoint=[0, 0, 0]
spacingX=1
spacingY=1
spacingZ=1

In this way, x=i, y=j, z=k. Since this image is a 2D image, the k should be 0 and 'z' should also be 0.

Then, I rotate the image, z is not 0. Like the following picture.

I don't know why z is -0.24.

It means the following conversion is wrong. And how can I obtain the image index by the world point?

worldPoint (x,y,z) = i*spacingX*directionX + j*spacingY*directionY + k*spacingZ*directionZ + originPoint

Any suggestion is appreciated!

Answer 1

vtkImageData has the method TransformPhysicalPointToContinuousIndex for going from world space to image space and TransformIndexToPhysicalPoint to go the other way.

I don't think the computation you're doing is right, since direction is 3x3 rotation matrix.