时间戳

首页 > 解决方案 > 如何在python中获得dicom切片的3D重建？

3d - 如何在python中获得dicom切片的3D重建？

问题描述

我有一个包含单次扫描的.dcm文件的目录。我能够获得 2D 视图（我与其他 dicom 查看器进行了交叉检查）。但我无法让 3D 视图正常工作。我尝试使用vtkDICOMImageReader该类，但它无法读取文件。所以我尝试从 3D numpy 数组中获取一个 Volume 对象以使用vtkplotter. 出来的观点显然是错误的。我认为 3D 数组需要一些处理。

import time
import glob
import pydicom
import numpy as np
from vtkplotter import Volume
import sys, os

def main(folderPath):
    st = time.time()
    my_glob = glob.glob1(folderPath, "*")
    numFiles = 0
    rejected = 0

    # return if empty directory
    if len(my_glob) == 0:
        return False

    # get all readable dicom files in  array
    tem = []
    for file in list(my_glob):
        try:
            data_item = pydicom.dcmread(os.path.join(folderPath, file))
            if hasattr(data_item, 'SliceLocation'):
                tem.append(data_item)
                numFiles += 1
            else:
                rejected += 1
                print(file)
        except Exception as e:
            pass
    print("read done %s | %d files | %d rejected" % (time.time() - st, numFiles, rejected))
    if len(tem) <= 0:
        return False

    tem.sort(key=lambda x: x.InstanceNumber)

    # make 3d np array from all slices
    unset = True
    for i in range(len(tem)):
        arr = tem[i].pixel_array.astype(np.float32)
        if unset:
            imShape = (arr.shape[0], arr.shape[1], len(tem))
            scaledIm = np.zeros(imShape)
            pix_spacing = tem[i].PixelSpacing
            dist = 0
            for j in range(2):
                cs = [float(q) for q in tem[j].ImageOrientationPatient]
                ipp = [float(q) for q in tem[j].ImagePositionPatient]
                parity = pow(-1, j)
                dist += parity*(cs[1]*cs[5] - cs[2]*cs[4])*ipp[0]
                dist += parity*(cs[2]*cs[3] - cs[0]*cs[5])*ipp[1]
                dist += parity*(cs[0]*cs[4] - cs[1]*cs[3])*ipp[2]
            z_spacing = abs(dist)
            slope = tem[i].RescaleSlope
            intercept = tem[i].RescaleIntercept
            unset = False
        scaledIm[:, :, i] = arr

    # convert to hounsfield units
    scaledIm = slope*scaledIm + intercept
    pix_spacing.append(z_spacing)

    wl = 300    # window parameters for Angio
    ww = 600

    windowed = np.zeros(imShape, dtype=np.uint8)
    # allImages[scaledIm <= (wl-0.5-(ww-1)/2.0)] = 0
    k = np.logical_and(scaledIm > (wl-0.5-(ww-1)/2.0), scaledIm <= (wl-0.5+(ww-1)/2.0))
    windowed[k] = ((scaledIm[k] - (wl-0.5))/(ww-1)+0.5)*255
    windowed[scaledIm > (wl-0.5+(ww-1)/2.0)] = 255
    # windowed image (in 2D) is correct i checked visually in other DICOM viewers
    print("arrays made %s" % (time.time() - st))


    # Volume(scaledIm, spacing=pix_spacing).show(bg="black")
    Volume(windowed, spacing=pix_spacing).show(bg="black")

    # X, Y, Z = np.mgrid[:30, :30, :30]
    # scalar_field = ((X-15)**2 + (Y-15)**2 + (Z-15)**2)/225
    # Volume(scalar_field, spacing=pix_spacing).show(bg="black")      # looks good on this example


if __name__ == '__main__':
    folder = sys.argv[1]
    main(folder)

需要做什么才能获得其他 dicom 查看器中显示的正确 3D 视图？

标签： 3dvtkdicompydicom

解决方案

我无法重现该问题，因为我收到来自 pydicom 的错误消息：

AttributeError: 'FileDataset' object has no attribute 'RescaleSlope'

无论如何，您可以尝试以下方法：

  • 更新到最新提交pip install git+https://github.com/marcomusy/vtkplotter.git

  • Volume将您的实例化修改为：

    # NOTE that pix_spacing[0] and pix_spacing[2] might be inverted
    vol = Volume(windowed, spacing=pix_spacing)
    vol.permuteAxes(2,1,0).mirror("y")
    vol.show(bg="black")

你也可以看看这个例子。

推荐阅读