使用python将dicom格式的文件转换为点云PLY

有时候我们会需要将CBCT扫描出来的dicom序列文件转换为点云格式，并进行下一步的处理。
那么这里我提供了一种思路，将dicom数据转换成了点云，但是这种转换仅保留了点的位置信息。可用于进一步的特征提取或者配准等工作

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)

结果展示：