Применить вращение, заданное углами Эйлера, к трехмерному изображению, в python - PullRequest
Новая
       79

Применить вращение, заданное углами Эйлера, к трехмерному изображению, в python

0 голосов
/ 14 января 2020

Я работаю с 3D-изображениями и должен поворачивать их в соответствии с углами Эйлера (phi, psi, theta) в соглашении 'zxz' (эти углы Эйлера являются частью набора данных, поэтому я должен использовать это соглашение). Я нашел функцию scipy.ndimage.rotate, которая кажется полезной в этом отношении. 

arrayR = scipy.ndimage.rotate(array , phi, axes=(0,1), reshape=False)
arrayR = scipy.ndimage.rotate(arrayR, psi, axes=(1,2), reshape=False)
arrayR = scipy.ndimage.rotate(arrayR, the, axes=(0,1), reshape=False)

К сожалению, это не делает то, что предполагалось. Вот почему:

Определение:

В соглашении zxz кадр xyz поворачивается три раза: сначала вокруг оси z на угол phi; затем о новой оси x на угол psi; затем о новейшей оси z на угол тета.

Однако с использованием приведенного выше кода повороты всегда выполняются относительно исходных осей. Вот почему полученные вращения не верны. У кого-нибудь есть предложение получить правильные повороты, как объяснено в определении?

Другими словами, в настоящем соглашении zxz повороты являются внутренними c (повороты вокруг осей вращающейся системы координат XYZ , солидарная с движущимся телом, которое меняет свою ориентацию после каждого элементарного вращения). Если я использую приведенный выше код, повороты будут extrinsi c (повороты вокруг осей xyz исходной системы координат, которая, как предполагается, остается неподвижной). Мне нужен способ для вращения extrinsi c, в python.

1 Ответ

0 голосов
/ 20 января 2020

Я нашел удовлетворительное решение по этой ссылке: https://nbviewer.jupyter.org/gist/lhk/f05ee20b5a826e4c8b9bb3e528348688

Этот метод использует np.meshgrid, scipy.ndimage.map_coordinates. Приведенная выше ссылка использует некоторую стороннюю библиотеку для генерации матрицы вращения, однако я использую scipy.spatial.transform.Rotation. Эта функция позволяет определять как внутренние вращения c, так и внешние вращения c: см. Описание scipy.spatial.transform.Rotation.from_euler.

Вот моя функция: 

import numpy as np
from scipy.spatial.transform import Rotation as R
from scipy.ndimage import map_coordinates

# Rotates 3D image around image center
# INPUTS
#   array: 3D numpy array
#   orient: list of Euler angles (phi,psi,the)
# OUTPUT
#   arrayR: rotated 3D numpy array
# by E. Moebel, 2020
def rotate_array(array, orient):
    phi = orient[0]
    psi = orient[1]
    the = orient[2]

    # create meshgrid
    dim = array.shape
    ax = np.arange(dim[0])
    ay = np.arange(dim[1])
    az = np.arange(dim[2])
    coords = np.meshgrid(ax, ay, az)

    # stack the meshgrid to position vectors, center them around 0 by substracting dim/2
    xyz = np.vstack([coords[0].reshape(-1) - float(dim[0]) / 2,  # x coordinate, centered
                     coords[1].reshape(-1) - float(dim[1]) / 2,  # y coordinate, centered
                     coords[2].reshape(-1) - float(dim[2]) / 2])  # z coordinate, centered

    # create transformation matrix
    r = R.from_euler('zxz', [phi, psi, the], degrees=True)
    mat = r.as_matrix()

    # apply transformation
    transformed_xyz = np.dot(mat, xyz)

    # extract coordinates
    x = transformed_xyz[0, :] + float(dim[0]) / 2
    y = transformed_xyz[1, :] + float(dim[1]) / 2
    z = transformed_xyz[2, :] + float(dim[2]) / 2

    x = x.reshape((dim[1],dim[0],dim[2]))
    y = y.reshape((dim[1],dim[0],dim[2]))
    z = z.reshape((dim[1],dim[0],dim[2])) # reason for strange ordering: see next line

    # the coordinate system seems to be strange, it has to be ordered like this
    new_xyz = [y, x, z]

    # sample
    arrayR = map_coordinates(array, new_xyz, order=1)

Примечание: Вы также можете использовать эту функцию для внутренних вращений c, просто приспособив первый аргумент 'from_euler' к вашему соглашению Эйлера. В этом случае вы получите эквивалентный результат, чем в моем первом посте (используя scipy.ndimage.rotate). Однако я заметил, что текущий код в 3 раза быстрее (0,01 с для тома 40 ^ 3), чем при использовании scipy.ndimage.rotate (0,03 с для тома 40 ^ 3).

Надеюсь, это кому-нибудь поможет!

...