CBCT是Cone Beam CT的简称,中文称为锥形束CT。它是一种口腔影像技术,利用锥形束投照计算机重组断层影像设备来获取三维图像。其原理是X线发生器以较低的射线量围绕投照体做环形DR(数字式投照),然后将多次数字投照后“交集”中所获得的数据在计算机中“重组”,进而获得三维图像。
要使用Python打开一个CBCT的口腔影像,你通常需要一个能够处理DICOM格式文件的库,因为DICOM是医学影像中常用的标准格式。pydicom和SimpleITK是两个常用的Python库,可以用来读取和处理DICOM文件。
首先,我们需要将这两个库安装一下:
pip install pydicom
pip install SimpleITK
1、查看CBCT切片
以下是一个使用pydicom库来读取DICOM文件的简单示例:
import pydicom
# 指定DICOM文件的路径
dicom_path = 'path_to_your_cbct_image.dcm'
# 读取DICOM文件
ds = pydicom.dcmread(dicom_path)
# 显示DICOM文件的信息
print(ds)
# 如果要查看图像数据,可以使用numpy库
import numpy as np
# 将DICOM数据转换为numpy数组
image_array = ds.pixel_array
# 你可以使用matplotlib等库来显示图像
import matplotlib.pyplot as plt
plt.imshow(image_array, cmap=plt.cm.bone) # 使用适合医学影像的colormap
plt.show()
这样我们就能够看到一个CT数据的切片了。但是如果我们打印ds里面的数据我们就会发现dicom数据范围往往不是0-255之间的数据。
由于医学成像设备(如CT、MRI等)的工作原理和人体组织的不同特性,医学图像的像素值(或称为灰度值)可以覆盖一个相对较大的范围。例如,空气的像素值可能接近于-1000,而骨头的像素值可能高达1000左右。这种范围的变化反映了不同组织对X射线或其他成像技术的吸收和散射特性的差异。
所以我们通常会对dicom数据进行一些预处理:
import matplotlib.pyplot as plt
import pydicom
import numpy as np
# 读取DICOM文件
dicom_path = 'path_to_your_dicom_image.dcm'
ds = pydicom.dcmread(dicom_path)
# 获取像素数据
pixel_data = ds.pixel_array
# 获取窗宽和窗位(如果DICOM头中有这些信息)
window_width = ds.WindowWidth if 'WindowWidth' in ds else 1500
window_center = ds.WindowCenter if 'WindowCenter' in ds else 50
# 窗宽窗位调整
min_value = window_center - window_width // 2
max_value = window_center + window_width // 2
scaled_pixel_data = np.clip(pixel_data, min_value, max_value)
scaled_pixel_data = (scaled_pixel_data - min_value) / (max_value - min_value) * 255
scaled_pixel_data = scaled_pixel_data.astype(np.uint8)
# 显示图像
plt.imshow(scaled_pixel_data, cmap=plt.cm.bone)
plt.axis('off') # 不显示坐标轴
plt.show()
2、查看一个DICOM的文件序列
如果你有一个包含多个DICOM文件的序列(例如,一个CBCT扫描通常包含多个层面的图像),你可能需要遍历这些文件并将它们组合成一个三维图像。这可以通过读取每个DICOM文件的元数据来实现,特别是ImagePositionPatient和ImageOrientationPatient,这些元数据可以帮助你正确地定位每个层面的图像。
SimpleITK库提供了更高级的功能,可以直接处理DICOM序列,并将其加载为三维图像:
import SimpleITK as sitk
# 指定DICOM序列的目录路径
dicom_dir = 'path_to_your_cbct_series'
# 使用SimpleITK读取DICOM序列
reader = sitk.ImageSeriesReader()
dicom_names = reader.GetGDCMSeriesFileNames(dicom_dir)
image3D = reader.Execute(dicom_names)
# 显示图像
# sitk.Show(image3D, SimpleITK.sitkInt16)
请确保将path_to_your_cbct_image.dcm和path_to_your_cbct_series替换为你的CBCT DICOM文件或DICOM序列的实际路径。
如果你遇到任何问题,比如文件路径错误、文件损坏或格式不正确等,Python会抛出异常。确保你的DICOM文件是有效的,并且你有足够的权限来读取它们。
3、将dicom数据转为点云
有时候我们需要把dicom的数据格式转换成点云数据,这里的点云转换是我自己的想法。读取每一个dicom序列,并将dicom切片中的数值转换成0-255的数值,然后对大于某个阈值的值的坐标进行写入。
import pydicom
import numpy as np
import os
from tqdm import tqdm
import time
# 这里是使用生成一个ASCII格式的PLY文件,包含了点的XYZ坐标以及可选的颜色信息。如果你的点云没有颜色信息,只需要调用create_ply_file(points)即可。
# 注意:对于大型点云数据集,通常会考虑使用二进制PLY格式以减小文件大小和提高读写效率。上述例子仅针对ASCII PLY格式,转换到二进制格式需要修改写入部分,并使用struct.pack()进行数据打包。
# 这里寻找点云的点
def get_pointlist_from_dicom(dicom_dir, thresh_hold=0):
print("Loading dicom series".center(120, "="))
point_list = np.empty([0,3])
for dcm_f in tqdm(os.listdir(dicom_dir)):
file_path = dicom_dir + dcm_f
dcm = pydicom.dcmread(file_path)
dcm_array = dcm.pixel_array
# 获取窗宽和窗位(如果DICOM头中有这些信息)
window_width = dcm.WindowWidth if 'WindowWidth' in dcm else 1500
window_center = dcm.WindowCenter if 'WindowCenter' in dcm else 50
series_number = float(dcm.InstanceNumber)
# 窗宽窗位调整
min_value = window_center - window_width // 2 # -1024
max_value = window_center + window_width // 2 # 3072
scaled_pixel_data = np.clip(dcm_array, min_value, max_value)
scaled_pixel_data = (scaled_pixel_data - min_value) / (max_value - min_value) * 255
scaled_pixel_data = scaled_pixel_data.astype(np.uint8)
# 使用numpy
x_index, y_index = np.where(scaled_pixel_data>thresh_hold)
z_index = np.array([series_number]*len(x_index))
layer_points = np.concatenate([[x_index], [y_index], [z_index]], axis=0).T
point_list = np.concatenate([point_list, layer_points], axis=0)
# # 垃圾写法 比上面的写法慢了 34 倍
# for x in range(len(scaled_pixel_data)):
# for y in range(len(scaled_pixel_data)):
# if scaled_pixel_data[x][y] > thresh_hold:
# point_list.append([x, y, layer_index])
print("\n")
return point_list
# 将点写成点云文件
def write_ply_file(point_list, point_cloud_file):
print("Writting PLY file".center(120, "="))
with open(point_cloud_file, "w+") as f:
# 写入PLY文件头
f.write("ply\n")
f.write("format ascii 1.0\n")
f.write("element vertex %d\n" % len(point_list))
f.write("property float x\n")
f.write("property float y\n")
f.write("property float z\n")
f.write("property uchar red\n")
f.write("property uchar green\n")
f.write("property uchar blue\n")
f.write("end_header\n")
for (x,y,z) in tqdm(point_list):
f.write("%f %f %f %d %d %d\n" % (x, y, z, 0, 255, 0))
print("Finished!\n")
if __name__ == "__main__":
dicom_dir = "./publick/"
point_cloud_file = f"./output/tooth_{str(int(time.time()))}.ply"
if not os.path.exists("./output/"):
os.mkdir("./output/")
point_list = get_pointlist_from_dicom(dicom_dir, thresh_hold=125)
write_ply_file(point_list,point_cloud_file)
结果展示:
网友评论