python打开读取dicom（CBCT）的数据

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)

结果展示：