标签：voxel slice patient python spacing dataset cosine np datasets

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 （还未实现，待补充）

标签：voxel,slice,patient,python,spacing,dataset,cosine,np,datasets
来源： https://blog.csdn.net/normol/article/details/90936225

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。