编程语言
首页 > 编程语言> > voxel 转换为 patient coordinate(python实现)

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 (还未实现,待补充)

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