voxel 转换为 patient coordinate（python实现）

dcm

其实自己感觉还未完全理解（博客内容若有错误请指出），先记下来，等答辩、课题等事情弄好再重新学习并补充。

一些基础概念别人博客已经写的很好了，我理解的关键点为：
1、病人坐标系的xyz定义方向为LPS(并非所有的，一些集成3D slice的软件用的是RAS）
2、图像坐标xy定义方向：(0,0)代表图像左上角，（max(x),max(y)）代表图像右下角（这儿x的延伸方向也是i的延伸方向，y的延伸方向等同于j的延伸方向）
3、仿射变换矩阵定义为

4、最需要注意的是，在下面_merge_slice_pixel_arrays函数里也有注释说明，python中数组默认order=C，即变换最快的轴向反而靠后，也就是我们一般读序列进来会是(Z,Y,X)，而matlab数组的默认order=F，读进来是正常的(X,Y,Z)。
对于dcm数据中的pixel_data而言，其为(height,width)，用数组的index取值时
pixel_data[i][j]代表的是第i行，第j列，这儿与第三点仿射变换矩阵中的i，j含义对调了，因此下面的代码里原作者将order改为了F，且将pixel_data进行了转置，但由于这种操作会不利于后续的操作。我认为只要清楚了这点，在后续操作需要用到变换矩阵时，手动进行矩阵的transpose即可。

import pydicom
import numpy as np
import os

"""
DICOM voxel to patient coordinate system mapping

https://github.com/innolitics/dicom-numpy/tree/master/dicom_numpy
"""
path='C:\\Users\\dell7920\\Desktop\\object\\3Dircadb1.1\\PATIENT_DICOM\\'
list_of_dicom_files = os.listdir(path)
datasets = [pydicom.dcmread(path + f) for f in list_of_dicom_files]

def _requires_rescaling(dataset):
    return hasattr(dataset, 'RescaleSlope') or hasattr(dataset, 'RescaleIntercept')

def combine_slices(slice_datasets, rescale=None):
    voxels = _merge_slice_pixel_arrays(slice_datasets, rescale)
    transform = _ijk_to_patient_xyz_transform_matrix(slice_datasets)

    return voxels, transform

def _merge_slice_pixel_arrays(slice_datasets, rescale=None):
    first_dataset = slice_datasets[0]
    num_rows = first_dataset.Rows
    num_columns = first_dataset.Columns
    num_slices = len(slice_datasets)

    sorted_slice_datasets = _sort_by_slice_spacing(slice_datasets)
    # print(sorted_slice_datasets)

    if rescale is None:
        rescale = any(_requires_rescaling(d) for d in sorted_slice_datasets)

    if rescale:
        # 注意，numpy默认的order是C，matlab才是fortan，Fortan最先读入增长最快的那个轴向，ct中是(X,Y,Z)
        voxels = np.empty((num_rows, num_columns, num_slices), dtype=np.float32)
        for k, dataset in enumerate(sorted_slice_datasets):
            slope = float(getattr(dataset, 'RescaleSlope', 1))
            intercept = float(getattr(dataset, 'RescaleIntercept', 0))
            # 原githubvoxel用的order='F', dataset.pixel_array.T
            voxels[:, :, k] = dataset.pixel_array.astype(np.float32) * slope + intercept
    else:
        dtype = first_dataset.pixel_array.dtype
        voxels = np.empty((num_rows, num_columns, num_slices), dtype=dtype)
        for k, dataset in enumerate(sorted_slice_datasets):
        # 原githubvoxel用的order='F', dataset.pixel_array.T
            voxels[:, :, k] = dataset.pixel_array

    return voxels

def _sort_by_slice_spacing(slice_datasets):
    slice_spacing = _slice_positions(slice_datasets)
    # print("slice_spacing ",sorted(slice_spacing))
    return [d for (s, d) in sorted(zip(slice_spacing, slice_datasets))]

def _extract_cosines(image_orientation):
    row_cosine = np.array(image_orientation[:3])
    column_cosine = np.array(image_orientation[3:])
    slice_cosine = np.cross(row_cosine, column_cosine)
    # print("slice cosine:",slice_cosine)
    return row_cosine, column_cosine, slice_cosine

def _slice_positions(slice_datasets):
    """
    获取切片在patient coordinate下的z轴坐标
    """
    image_orientation = slice_datasets[0].ImageOrientationPatient
    row_cosine, column_cosine, slice_cosine = _extract_cosines(image_orientation)
    return [np.dot(slice_cosine, d.ImagePositionPatient) for d in slice_datasets]


def _slice_spacing(slice_datasets):
    """
    获取切片在patient coordinate下的平均间隔
    """
    if len(slice_datasets) > 1:
        slice_positions = _slice_positions(slice_datasets)
        slice_positions_diffs = np.diff(sorted(slice_positions))
        return np.mean(slice_positions_diffs)
    else:
        return 0.0

def _ijk_to_patient_xyz_transform_matrix(slice_datasets):
    first_dataset = _sort_by_slice_spacing(slice_datasets)[0]
    image_orientation = first_dataset.ImageOrientationPatient
    row_cosine, column_cosine, slice_cosine = _extract_cosines(image_orientation)

    row_spacing, column_spacing = first_dataset.PixelSpacing
    slice_spacing = _slice_spacing(slice_datasets)

    transform = np.identity(4, dtype=np.float32)

    transform[:3, 0] = row_cosine * column_spacing
    transform[:3, 1] = column_cosine * row_spacing
    transform[:3, 2] = slice_cosine * slice_spacing

    transform[:3, 3] = first_dataset.ImagePositionPatient

    return transform

voxels, transform = combine_slices(datasets)

print(transform)

参考：
http://nipy.org/nipy/devel/code_discussions/understanding_affines.html
https://nben.net/MRI-Geometry/#mri-geometry
https://brainder.org/2012/09/23/the-nifti-file-format/
https://blog.csdn.net/zssureqh/article/details/61636150

nii （还未实现，待补充）