DICOM 文件读取与显示 - Python 实现
import cv2
import numpy as np
# 定义原始像素类型
class E_SourcePixelType:
ESourcePixelType_U16 = 0 # USHORT
ESourcePixelType_I16 = 1 # SHORT
# 读取Dicom文件
def read_dicom_file(file_path):
# 初始化变量
is_VR = True
is_LittleEndian = True
file_length = 0
VR = [0, 0, 0]
pixDataLen = 0
pixDataOffset = 0
channel = 0
rows = 0
cols = 0
dataLen = 0
validLen = 0
pixelType = None
windowsWidth = 0
windowsCenter = 0
ZeroIsBlack = True
RescaleSlope = 0.06
RescaleIntercept = 0
# 打开文件
f = open(file_path, "rb")
# 读取文件长度并返回文件起始位置
f.seek(0, 2)
file_length = f.tell()
f.seek(0, 0)
# 跳过128字节的DICOM文件头
f.seek(128, 0)
# 判断DICOM文件头是否为DICM
headchar = f.read(4)
if headchar != b"DICM":
f.close()
print("File is not DICM")
return None
# 读取DICOM文件头信息
while f.tell() + 6 < file_length:
# 读取Tag
tag = f.read(4)
len = 0
# 处理不同Tag类型
if tag[0] == 0x02:
VR = f.read(2)
if VR == b"OB" or VR == b"OW" or VR == b"SQ":
# 跳过2字节
f.seek(2, 1)
# 读取长度
len = int.from_bytes(f.read(4), byteorder="little")
else:
# 读取长度
len = int.from_bytes(f.read(2), byteorder="little")
elif tag[0] == 0xff and tag[1] == 0xfe:
if tag[2] == 0xe0 and tag[3] == 0x00:
# 读取长度
len = int.from_bytes(f.read(4), byteorder="little")
elif is_VR:
VR = f.read(2)
if VR == b"OB" or VR == b"OW" or VR == b"SQ":
# 跳过2字节
f.seek(2, 1)
# 读取长度
len = int.from_bytes(f.read(4), byteorder="little")
else:
# 读取长度
len = int.from_bytes(f.read(2), byteorder="little")
else:
# 读取长度
len = int.from_bytes(f.read(4), byteorder="little")
# 解析不同Tag的值
if tag[0] == 0x02 and tag[1] == 0x10:
msg = f.read(len)
# 判断字节序和VR
if msg == b"1.2.840.10008.1.2.1":
is_LittleEndian = True
is_VR = True
elif msg == b"1.2.840.10008.1.2.2":
is_LittleEndian = False
is_VR = True
elif msg == b"1.2.840.10008.1.2":
is_LittleEndian = True
is_VR = False
elif tag[0] == 0x28 and tag[1] == 0x103:
# 读取像素类型
m = int.from_bytes(f.read(2), byteorder="little")
if m == 0:
pixelType = E_SourcePixelType.ESourcePixelType_U16
elif m == 1:
pixelType = E_SourcePixelType.ESourcePixelType_I16
elif tag[0] == 0x7fe0 and tag[1] == 0x10:
# 读取像素数据长度和偏移量
pixDataLen = len
pixDataOffset = f.tell()
# 跳过像素数据
f.seek(len, 1)
elif tag[0] == 0x28 and tag[1] == 0x10:
# 读取行数
rows = int.from_bytes(f.read(2), byteorder="little")
elif tag[0] == 0x28 and tag[1] == 0x11:
# 读取列数
cols = int.from_bytes(f.read(2), byteorder="little")
elif tag[0] == 0x28 and tag[1] == 0x02:
# 读取通道数
channel = int.from_bytes(f.read(2), byteorder="little")
elif tag[0] == 0x28 and tag[1] == 0x101:
# 读取有效位数
validLen = int.from_bytes(f.read(2), byteorder="little")
elif tag[0] == 0x28 and tag[1] == 0x100:
# 读取数据长度
dataLen = int.from_bytes(f.read(2), byteorder="little")
elif tag[0] == 0x28 and tag[1] == 0x1050:
# 读取窗口中心值
msg = f.read(len)
windowsCenter = int(msg)
elif tag[0] == 0x28 and tag[1] == 0x1051:
# 读取窗口宽度值
msg = f.read(len)
windowsWidth = int(msg)
elif tag[0] == 0x0028 and tag[1] == 0x0004:
# 读取图像颜色类型
msg = f.read(len)
if msg == b"MONOCHROME1 ":
ZeroIsBlack = False
elif msg == b"MONOCHROME2 ":
ZeroIsBlack = True
elif tag[0] == 0x0028 and tag[1] == 0x1052:
# 读取重标定截距
msg = f.read(len)
RescaleIntercept = float(msg)
elif tag[0] == 0x0028 and tag[1] == 0x1053:
# 读取重标定斜率
msg = f.read(len)
RescaleSlope = float(msg)
else:
# 跳过其他Tag
f.seek(len, 1)
# 返回文件起始位置
f.seek(pixDataOffset, 0)
# 如果窗口中心和宽度为0，则计算默认值
if windowsCenter == 0 and windowsWidth == 0:
windowsWidth = 1 << validLen
windowsCenter = windowsWidth // 2
# 创建图像
src = np.zeros((rows, cols), dtype=np.uint8)
src2 = np.zeros((rows, cols), dtype=np.int16)
# 计算像素值转换参数
fCtA = 256 / windowsWidth
fCtB = 128 - 256 * windowsCenter / windowsWidth
if fCtB < 0:
fCtB = 0
if fCtB > 255:
fCtB = 255
# 处理不同通道数的图像
if channel == 1:
# 单通道图像
for i in range(rows):
for j in range(cols):
# 读取像素值
pix = f.read(2)
gray = 0
gray2 = 0
# 根据像素类型处理像素值
if pixelType == E_SourcePixelType.ESourcePixelType_U16:
if validLen <= 8:
if is_LittleEndian:
gray = pix[0]
else:
gray = pix[1]
else:
temp = 0
if is_LittleEndian:
gray = int.from_bytes(pix, byteorder="little")
if gray > 32767:
gray = 32767
temp = gray * RescaleSlope + RescaleIntercept
temp = temp * fCtA + fCtB
else:
gray = pix[1] + pix[0] * 256
temp = gray * RescaleSlope + RescaleIntercept
temp = temp * fCtA + fCtB
nValue = int(temp)
if nValue > 255:
nValue = 255
elif nValue < 0:
nValue = 0
nPixel = nValue
elif pixelType == E_SourcePixelType.ESourcePixelType_I16:
if validLen <= 8:
if is_LittleEndian:
gray2 = pix[0]
else:
gray2 = pix[1]
else:
temp = 0
if is_LittleEndian:
gray2 = int.from_bytes(pix, byteorder="little", signed=True)
if gray2 > 32767:
gray2 = 32767
if gray2 < -32767:
gray2 = -32767
temp = gray2 * RescaleSlope + RescaleIntercept
temp = temp * fCtA + fCtB
else:
gray2 = pix[1] + pix[0] * 256
temp = gray2 * RescaleSlope + RescaleIntercept
temp = temp * fCtA + fCtB
nValue = int(temp)
if nValue > 255:
nValue = 255
elif nValue < 0:
nValue = 0
nPixel = nValue
# 根据颜色类型调整像素值
if not ZeroIsBlack:
nPixel = 255 - nPixel
# 设置像素值
src[i, j] = nPixel
elif channel == 3:
# 三通道图像
for i in range(rows):
for j in range(cols):
# 读取像素值
pix = f.read(3)
# 设置像素值
src[i, j] = pix[::-1]
# 关闭文件
f.close()
# 显示图像
cv2.imshow("Dicomimage", src)
cv2.waitKey(0)
cv2.destroyAllWindows()
# 测试代码
file_path = "C:\\Users\\lenovo\\Desktop\\数据结构与算法C++\\20084125-张亭-数据结构算法读取Dicom文件\\11.dcm"
read_dicom_file(file_path)
原文地址: https://www.cveoy.top/t/topic/naBA 著作权归作者所有。请勿转载和采集!