使用pydicom实现Dicom文件读取与CT图像窗宽窗位调整

使用pydicom实现Dicom文件读取与CT图像窗宽窗位调整1.前言为了能够在Labelme上对Dicom图像进行编辑,这里对python环境下Dicom文件的读取进行了研究。在Dicom图像中CT的窗宽窗位是一个很重要的概念,但是找了半天在pydicom中没有相关设置函数,这里跟DCMTK还不一样。但是可以根据两个tag得到CT图像的CT值,那就是(0028|1052):rescaleintercept和(0028|1053):rescales…

大家好,又见面了,我是你们的朋友全栈君。

1. 前言

为了能够在Labelme上对Dicom图像进行编辑,这里对python环境下Dicom文件的读取进行了研究。在Dicom图像中CT的窗宽窗位是一个很重要的概念,但是找了半天在pydicom中没有相关设置函数,这里跟DCMTK还不一样。但是可以根据两个tag得到CT图像的CT值,那就是(0028|1052):rescale intercept和(0028|1053):rescale slope。则按照下面的算子得到CT图像,进而就可以调整窗宽窗位了

Hu = pixel * slope + intercept

至于那个部位的窗宽窗位是多少各位看官就可以自行百度了。

2. 代码实现

# -*- coding=utf-8 -*-
import matplotlib.pyplot as plt
import pydicom
import pydicom.uid
import sys
import PIL.Image as Image
from PyQt5 import QtGui
import os

step1:读取Dicom图像数据与得到CT值图像(CT图)

have_numpy = True

try:
    import numpy
except ImportError:
    have_numpy = False
    raise

sys_is_little_endian = (sys.byteorder == 'little')

NumpySupportedTransferSyntaxes = [
    pydicom.uid.ExplicitVRLittleEndian,
    pydicom.uid.ImplicitVRLittleEndian,
    pydicom.uid.DeflatedExplicitVRLittleEndian,
    pydicom.uid.ExplicitVRBigEndian,
]

# 支持的传输语法
def supports_transfer_syntax(dicom_dataset):
    """ Returns ------- bool True if this pixel data handler might support this transfer syntax. False to prevent any attempt to try to use this handler to decode the given transfer syntax """
    return (dicom_dataset.file_meta.TransferSyntaxUID in
            NumpySupportedTransferSyntaxes)


def needs_to_convert_to_RGB(dicom_dataset):
    return False


def should_change_PhotometricInterpretation_to_RGB(dicom_dataset):
    return False


# 加载Dicom图像数据
def get_pixeldata(dicom_dataset):
    """If NumPy is available, return an ndarray of the Pixel Data. Raises ------ TypeError If there is no Pixel Data or not a supported data type. ImportError If NumPy isn't found NotImplementedError if the transfer syntax is not supported AttributeError if the decoded amount of data does not match the expected amount Returns ------- numpy.ndarray The contents of the Pixel Data element (7FE0,0010) as an ndarray. """
    if (dicom_dataset.file_meta.TransferSyntaxUID not in
            NumpySupportedTransferSyntaxes):
        raise NotImplementedError("Pixel Data is compressed in a "
                                  "format pydicom does not yet handle. "
                                  "Cannot return array. Pydicom might "
                                  "be able to convert the pixel data "
                                  "using GDCM if it is installed.")

    # 设置窗宽窗位
    #dicom_dataset.

    if not have_numpy:
        msg = ("The Numpy package is required to use pixel_array, and "
               "numpy could not be imported.")
        raise ImportError(msg)
    if 'PixelData' not in dicom_dataset:
        raise TypeError("No pixel data found in this dataset.")

    # Make NumPy format code, e.g. "uint16", "int32" etc
    # from two pieces of info:
    # dicom_dataset.PixelRepresentation -- 0 for unsigned, 1 for signed;
    # dicom_dataset.BitsAllocated -- 8, 16, or 32
    if dicom_dataset.BitsAllocated == 1:
        # single bits are used for representation of binary data
        format_str = 'uint8'
    elif dicom_dataset.PixelRepresentation == 0:
        format_str = 'uint{}'.format(dicom_dataset.BitsAllocated)
    elif dicom_dataset.PixelRepresentation == 1:
        format_str = 'int{}'.format(dicom_dataset.BitsAllocated)
    else:
        format_str = 'bad_pixel_representation'
    try:
        numpy_dtype = numpy.dtype(format_str)
    except TypeError:
        msg = ("Data type not understood by NumPy: "
               "format='{}', PixelRepresentation={}, "
               "BitsAllocated={}".format(
                   format_str,
                   dicom_dataset.PixelRepresentation,
                   dicom_dataset.BitsAllocated))
        raise TypeError(msg)

    if dicom_dataset.is_little_endian != sys_is_little_endian:
        numpy_dtype = numpy_dtype.newbyteorder('S')

    pixel_bytearray = dicom_dataset.PixelData

    if dicom_dataset.BitsAllocated == 1:
        # if single bits are used for binary representation, a uint8 array
        # has to be converted to a binary-valued array (that is 8 times bigger)
        try:
            pixel_array = numpy.unpackbits(
                numpy.frombuffer(pixel_bytearray, dtype='uint8'))
        except NotImplementedError:
            # PyPy2 does not implement numpy.unpackbits
            raise NotImplementedError(
                'Cannot handle BitsAllocated == 1 on this platform')
    else:
        pixel_array = numpy.frombuffer(pixel_bytearray, dtype=numpy_dtype)
    length_of_pixel_array = pixel_array.nbytes
    expected_length = dicom_dataset.Rows * dicom_dataset.Columns
    if ('NumberOfFrames' in dicom_dataset and
            dicom_dataset.NumberOfFrames > 1):
        expected_length *= dicom_dataset.NumberOfFrames
    if ('SamplesPerPixel' in dicom_dataset and
            dicom_dataset.SamplesPerPixel > 1):
        expected_length *= dicom_dataset.SamplesPerPixel
    if dicom_dataset.BitsAllocated > 8:
        expected_length *= (dicom_dataset.BitsAllocated // 8)
    padded_length = expected_length
    if expected_length & 1:
        padded_length += 1
    if length_of_pixel_array != padded_length:
        raise AttributeError(
            "Amount of pixel data %d does not "
            "match the expected data %d" %
            (length_of_pixel_array, padded_length))
    if expected_length != padded_length:
        pixel_array = pixel_array[:expected_length]
    if should_change_PhotometricInterpretation_to_RGB(dicom_dataset):
        dicom_dataset.PhotometricInterpretation = "RGB"
    if dicom_dataset.Modality.lower().find('ct') >= 0:  # CT图像需要得到其CT值图像
        pixel_array = pixel_array * dicom_dataset.RescaleSlope + dicom_dataset.RescaleIntercept  # 获得图像的CT值
    pixel_array = pixel_array.reshape(dicom_dataset.Rows, dicom_dataset.Columns*dicom_dataset.SamplesPerPixel)
    return pixel_array, dicom_dataset.Rows, dicom_dataset.Columns

step2:对于CT图像设置窗宽窗位

# 调整CT图像的窗宽窗位
def setDicomWinWidthWinCenter(img_data, winwidth, wincenter, rows, cols):
    img_temp = img_data
    img_temp.flags.writeable = True
    min = (2 * wincenter - winwidth) / 2.0 + 0.5
    max = (2 * wincenter + winwidth) / 2.0 + 0.5
    dFactor = 255.0 / (max - min)

    for i in numpy.arange(rows):
        for j in numpy.arange(cols):
            img_temp[i, j] = int((img_temp[i, j]-min)*dFactor)

    min_index = img_temp < 0
    img_temp[min_index] = 0
    max_index = img_temp > 255
    img_temp[max_index] = 255

    return img_temp

step3:获取Dicom中的tag信息

第一种方式:

# 加载Dicom图片中的Tag信息
def loadFileInformation(filename):
    information = {}
    ds = pydicom.read_file(filename)
    information['PatientID'] = ds.PatientID
    information['PatientName'] = ds.PatientName
    information['PatientBirthDate'] = ds.PatientBirthDate
    information['PatientSex'] = ds.PatientSex
    information['StudyID'] = ds.StudyID
    information['StudyDate'] = ds.StudyDate
    information['StudyTime'] = ds.StudyTime
    information['InstitutionName'] = ds.InstitutionName
    information['Manufacturer'] = ds.Manufacturer
    print(dir(ds))
    print(type(information))
    return information

第二种方式

dcm = pydicom.dcmread(fileanme)  # 加载Dicom数据

print(dcm[0x0008, 0x0060])
>>(0008, 0060) Modality                            CS: 'MR'
print(dcm[0x0008, 0x0060].VR)
>>CS
print(dcm[0x0008, 0x0060].value)
>>MR

step4:Dicom图像数据转换为PIL.Image

dcm = pydicom.dcmread(fileanme)  # 加载Dicom数据
dcm_img = Image.fromarray(img_data)  # 将Numpy转换为PIL.Image
dcm_img = dcm_img.convert('L')

# 保存为jpg文件,用作后面的生成label用
dcm_img.save('temp.jpg')
# 显示图像
dcm_img.show()

3. 结果展示

调整了窗宽窗位的脑部CT图像:
这里写图片描述

4. 参考资料

  1. Pydicom User Guide
  2. 【医学影像】窗宽窗位与其处理方法
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请联系我们举报,一经查实,本站将立刻删除。

发布者:全栈程序员-站长,转载请注明出处:https://javaforall.net/148747.html原文链接:https://javaforall.net

(1)
全栈程序员-站长的头像全栈程序员-站长


相关推荐

  • Java 排列组合_c语言排列组合函数

    Java 排列组合_c语言排列组合函数importjava.util.Arrays;//利用二进制算法进行全排列//count1:170187//count2:291656publicclasstest{publicstaticvoidmain(String[]args){longstart=System.currentTimeMillis();count2();longend=System.currentTim…

    2022年8月30日
    2
  • Teahour 以太坊专访-文字版

    Teahour 以太坊专访-文字版原文地址:http://ethfans.org/shaoping/articles/talk-with-jan-about-ehtereum本文是播客Teahour以太坊专访文字版,收听地址 http://teahour.fm/2016/01/19/talk-with-jan-about-ehtereum.html ,也可以在喜马拉雅、iTunes搜索Teahour收听。(前边的开

    2022年5月16日
    40
  • 移动apn接入点哪个快(移动哪个接入点网速快)

    9条解答1.中国移动4g接入点最快我不同意,中国LTE才是最快的,就是流量怕你不够,名称:LTE,APN:lte(小写的),APN协议漫游协议lpv4,我的达到了2.5~3M/S,骗你出家2.4g网速很慢,求助apn设置方法您好!您先进入手机的接入点设置–新建apn接入点–名称乱填,apn:cmtds–保存保存之后,选择自己刚刚设置的接入点,然后您就会发现您自己的手机4G网络的速度很快了,…

    2022年4月12日
    465
  • XSHELL安装指南

    XSHELL安装指南开发环境部署目的:利用ssh远程登陆服务器(在windows系统下远程连接linux)下载XSHELL7XSHELL7下载网址:https://www.netsarang.com/zh/xshell/点击“下载”点击“免费授权界面”以上是XSHELL7的下载过程然后找到右键“以管理员身份运行”一上来会出现这种错误,先点击“是(Y)”过程中一直点击“下一步”,以及“我同意”类似的,然后选择个安装路径就可以没啥特殊的。到最后一切顺利的话会显示下面这样的界面一般通向成功的道

    2025年6月17日
    0
  • tensorflow2.0手写数字识别_python 数字识别

    tensorflow2.0手写数字识别_python 数字识别本文使用Tensorflow框架进行Python编程实现基于卷积神经网络的手写数字识别算法,并将其封装在一个GUI界面中,最终,设计并实现了一个手写数字识别系统。

    2022年9月14日
    2
  • 回归分析详解及matlab实现[通俗易懂]

    回归分析详解及matlab实现[通俗易懂] 回归分析方法 想要资源的请关注公众号:在一起的足球自动获取资源和数十种经典算法,帮助各位提升自己之前留的是自己的qq号感觉好多资源都不能让大家自行选择本着开源的精神,在公众号挂了百度云链接,这样晚上就不用被吵醒了当人们对研究对象的内在特性和各因素间的关系有比较充分的认识时,一般用机理分析方法建立数学模型。如果由于客观事物内部规律的复杂性及人们认识程度的限制,无法分析实际对象内在的因果关系,…

    2025年7月20日
    5

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

关注全栈程序员社区公众号