如何匹配 vtk polydata 和 itk 转换

Question

我需要使用相同的变换矩阵变换（暂时旋转）itk 图像和 vtk 多数据，但我遇到了麻烦。

所有代码和测试数据都在这里：https://github.com/jmerkow/vtk_itk_rotate

相关部分如下：

import SimpleITK as sitk
import vtk
import numpy as np
def rotate_img(img, rotation_center=None, theta_x=0,theta_y=0, theta_z=0, translation=(0,0,0), interp=sitk.sitkLinear, pixel_type=None, default_value=None):
    if not rotation_center:
        rotation_center = np.array(img.GetOrigin())+np.array(img.GetSpacing())*np.array(img.GetSize())/2
    if default_value is None:
        default_value = img.GetPixel(0,0,0)
    pixel_type = img.GetPixelIDValue()

    rigid_euler = sitk.Euler3DTransform(rotation_center, theta_x, theta_y, theta_z, translation)
    return sitk.Resample(img, img, rigid_euler, interp, default_value, pixel_type)

def rotate_polydata(pd, rotation_center, theta_x=0,theta_y=0, theta_z=0, translation=(0,0,0)):
    rigid_euler = sitk.Euler3DTransform(rotation_center, -theta_x, -theta_y, -theta_z, translation)
    matrix = np.zeros([4,4])
    old_matrix=np.array(rigid_euler.GetMatrix()).reshape(3,3)
    matrix[:3,:3] = old_matrix
    matrix[-1,-1] = 1

    # to rotate about a center we first need to translate
    transform_t = vtk.vtkTransform()
    transform_t.Translate(-rotation_center)
    transformer_t = vtk.vtkTransformPolyDataFilter()
    transformer_t.SetTransform(transform_t)
    transformer_t.SetInputData(pd)
    transformer_t.Update()

    transform = vtk.vtkTransform()
    transform.SetMatrix(matrix.ravel())

    transformer = vtk.vtkTransformPolyDataFilter()
    transformer.SetTransform(transform)
    transformer.SetInputConnection(transformer_t.GetOutputPort())
    transformer.Update()

    # translate back
    transform_t2 = vtk.vtkTransform()
    transform_t2.Translate(rotation_center)
    transformer_t2 = vtk.vtkTransformPolyDataFilter()
    transformer_t2.SetTransform(transform_t2)
    transformer_t2.SetInputConnection(transformer.GetOutputPort())
    transformer_t2.Update()

    return transformer_t2.GetOutputDataObject(0)

datafn = 'test.mha'
polydata_file = 'test.vtp'
reader = vtk.vtkXMLPolyDataReader()
reader.SetFileName(polydata_file)
reader.Update()
pd = reader.GetOutput()

img = sitk.ReadImage(datafn)
seg = pd_to_itk_image(pd, img)
rotation_center = np.array(img.GetOrigin())+np.array(img.GetSpacing())*np.array(img.GetSize())/2
thetas = [0, 50]
thetas = [0, 50]
for theta_x in thetas:
    for theta_y in thetas:
        for theta_z in thetas:
            theta_xr = theta_x/180.*np.pi
            theta_yr = theta_y/180.*np.pi
            theta_zr = theta_z/180.*np.pi
            img_rot=rotate_img(img, theta_z=theta_zr, theta_y=theta_yr, theta_x=theta_xr)
            seg_rot=rotate_img(seg, theta_z=theta_zr, theta_y=theta_yr, theta_x=theta_xr, interp=sitk.sitkNearestNeighbor, default_value=0)
            pd_rot = rotate_polydata(pd, rotation_center, theta_z=theta_zr, theta_y=theta_yr, theta_x=theta_xr)
            seg_pd_rot = pd_to_itk_image(pd_rot, img_rot)
            mse = ((sitk.GetArrayFromImage(seg_pd_rot)-sitk.GetArrayFromImage(seg_rot))**2.).mean()

            print theta_x, theta_y, theta_z, mse

#this outputs for this particular volume:
#0 0 0 mse: 0.0
#0 0 50 mse: 50.133369863 visually about the same
#0 50 0 mse: 25.2197787166 visually about the same
#0 50 50 mse: 863.588476181 visually totally different
#50 0 0 mse: 20.4021692276 visually about the same
#50 0 50 mse: 546.699844301 visually totally different
#50 50 0 mse: 662.337975204 visually totally different
#50 50 50 mse: 339.220945537 visually totally different

此代码旋转从多数据生成的二进制体积，并对多数据执行相同的旋转操作，然后从中生成二进制体积。我希望这两个结果（大致）相同，但是，如果我围绕一个以上的轴旋转，我得到的是两个完全不同的旋转。这让我感到困惑，因为我从一个矩阵中获取变换矩阵并将其直接应用于另一个矩阵。

如何设置这些转换，使两个操作执行相同的转换？为什么我们最终会得到不同的结果？

Answer 1

欧拉角的顺序对最终结果很重要 [Wikipedia]. Also, matrix pre-multiplying has reversed order to post-multiplying [vtkTransform]。尝试调用 vtkTransform::PostMultiply() 或反转 rotate_polydata 函数中的转换顺序。这个很容易尝试。

如果这不能解决问题，请查看 ITK 如何在 ComputeOffset, TransformPoint and ComputeMatrix, and how VTK does it in vtkLinearTransformPoint 中应用转换。这应该可以解释行为上的差异并提供有关如何实现相同转换的线索。

Answer 2

感谢 Dženan 为我指明了正确的方向。

在这种情况下，答案很简单。 VTK 和 ITK 使用不同的 row/column 主要格式进行矩阵乘法。所以答案就是在将矩阵放入 vtkTransform 之前转置矩阵。

这是新功能。

def rotate_polydata(pd, rotation_center, theta_x=0,theta_y=0, theta_z=0):
    #I don't want to deal with translation
    translation=(0,0,0)
    rigid_euler = sitk.Euler3DTransform(rotation_center, theta_x, theta_y, theta_z, translation)
    matrix = np.zeros([4,4])
    old_matrix=np.array(rigid_euler.GetMatrix()).reshape(3,3)
    matrix[:3,:3] = old_matrix
    matrix[-1,-1] = 1
    #ITK and VTK use different orders.
    matrix= matrix.T

    # to rotate about a center we first need to translate
    transform_t = vtk.vtkTransform()
    transform_t.Translate(-rotation_center)
    transformer_t = vtk.vtkTransformPolyDataFilter()
    transformer_t.SetTransform(transform_t)
    transformer_t.SetInputData(pd)
    transformer_t.Update()

    transform = vtk.vtkTransform()
    transform.SetMatrix(matrix.ravel())
    transform.Translate(translation)
    transform.PostMultiply()

    transformer = vtk.vtkTransformPolyDataFilter()
    transformer.SetTransform(transform)
    transformer.SetInputConnection(transformer_t.GetOutputPort())
    transformer.Update()

    # translate back
    transform_t2 = vtk.vtkTransform()
    transform_t2.Translate(rotation_center)
    transformer_t2 = vtk.vtkTransformPolyDataFilter()
    transformer_t2.SetTransform(transform_t2)
    transformer_t2.SetInputConnection(transformer.GetOutputPort())
    transformer_t2.Update()

    return transformer_t2.GetOutputDataObject(0)

如何匹配 vtk polydata 和 itk 转换

how to match vtk polydata and itk transforms

python

vtk

itk