как совместить vtk polydata и itk transforms

Question

Мне нужно преобразовать (вращение на данный момент) изображение itk и vtk polydata с использованием той же матрицы преобразования, но у меня возникают проблемы.

Весь код и данные испытаний здесь: https://github.com/jmerkow/vtk_itk_rotate

Вот relavent части:

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 

Этот код поворачивает двоичный объем генерируемого из polydata, и выполняет ту же операцию вращения на polydata затем генерирует из него двоичный том. Я ожидаю, что эти два результата будут (примерно) одинаковыми, однако я получаю два совершенно разных вращения, если я вращаюсь вокруг более чем одной оси. Это меня озадачивает, так как я принимаю матрицу преобразования от одного и применяю ее непосредственно к другому.

Как настроить эти преобразования таким образом, чтобы обе операции выполняли одни и те же преобразования? И почему мы заканчиваем разными результатами?

Answer 1

Благодарим вас за Dženan за то, что указали мне нужное направление.

Ответ, в данном случае, был прост. VTK и ITK используют различные форматы строк/столбцов для их матричного умножения. Таким образом, ответ заключался просто в том, чтобы транспонировать матрицу, прежде чем вставлять ее в 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) 

Answer 2

Порядок углов Эйлера имеет значение для окончательного результата [Wikipedia]. Кроме того, предварительное умножение матрицы имеет обратный порядок для последующего умножения [vtkTransform]. Попробуйте позвонить vtkTransform::PostMultiply() или изменить порядок преобразований в функции rotate_polydata. Это легко попробовать.

Если это не решить, как смотреть в ИТК является применение преобразований в ComputeOffset, TransformPoint и ComputeMatrix, и как ВТК делает это в vtkLinearTransformPoint. Это должно объяснить различия в поведении и дать представление о том, как достичь тех же преобразований.