如何在python中获得dicom slice的3D重建?

Shrutesh Patil

我有一个包含单次扫描的.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视图,需要做什么?

us

我无法重现该问题,因为我从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")

您也可以查看此示例

本文收集自互联网,转载请注明来源。

如有侵权,请联系[email protected] 删除。

编辑于
0

我来说两句

0条评论
登录后参与评论

相关文章

来自分类Dev

如何在数字重建射线照相术中获得3D数据的一半

来自分类Dev

为R中的3D图重建数据框

来自分类Dev

如何在ThreeJS中获得3D对象尺寸?

来自分类Dev

使用Matlab中的内置功能进行3D重建(radon和iradon)

来自分类Dev

如何在C ++中的3D数组中获得三维数组?

来自分类Dev

如何在python中旋转3D文本?

来自分类Dev

如何在python中顺序填充3d矩阵?

来自分类Dev

如何在python中绘制3d图形

来自分类Dev

如何在Dlib C ++中获得3D头部姿势估计坐标轴

来自分类Dev

如何在通过kmeans获得的星团的R中绘制3D图?

来自分类Dev

3D重建和对象移动

来自分类Dev

3D重建和距离测量

来自分类Dev

如何在python中的2D列表中获得一行的平均值?

来自分类Dev

如何在将Matplotlib图形嵌入到Tkinter画布中时获得(3d)交互性

来自分类Dev

如何在python / plotly中制作2D向量分布的3D直方图

来自分类Dev

如何在Python中获得到加密数据库的SQLite3连接

来自分类Dev

如何在Python 3中获得组合Unicode字符的显示宽度?

来自分类Dev

如何在Python 3中获得以纳秒为单位的时间戳?

来自分类Dev

如何在sqlite3 python中获得更多的一个字母查询?

来自分类Dev

如何在python中转置3D列表?

来自分类Dev

如何在python的csv中导入3D表?

来自分类Dev

如何在Python上制作3D图形动画

来自分类Dev

如何在Unity 3d的场景中显示相机?

来自分类Dev

如何在Matplotlib中更新3D箭头动画

来自分类Dev

如何在Matlab中输入3D矩阵?

来自分类Dev

如何在R中绘制3D函数?

来自分类Dev

如何在Matlab中快速绘制3D向量?

来自分类Dev

如何在WPF 3d中调整模型大小?

来自分类Dev

如何在R中制作3D直方图

Related 相关文章

热门标签

归档