{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":71549,"databundleVersionId":8561470,"sourceType":"competition"}],"dockerImageVersionId":30698,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"一、导入库","metadata":{}},{"cell_type":"code","source":"# 从skimage.measure模块导入marching_cubes函数，用于从三维标量场中提取等值面  \nfrom skimage.measure import marching_cubes\n# 导入NumPy库，NumPy是Python的一个核心库，支持大量的维度数组与矩阵运算，  \n# 此外也提供了大量的数学函数库，是数据科学、科学计算和机器学习等领域的基础。  \nimport numpy as np  \n# 导入Pandas库，Pandas是一个强大的Python数据分析工具库，  \nimport pandas as pd  \n# 导入os库，os模块提供了许多与操作系统交互的功能，  \n# 比如文件和目录管理、环境变量访问、进程管理等。  \nimport os  \n# 导入time库，time模块提供了各种与时间相关的函数，  \n# 允许你获取当前时间、执行时间测量（如计算代码运行时间）等。  \nimport time  \n# 导入matplotlib.pyplot库，这是一个Python的绘图库，  \n# 它提供了一个类似于MATLAB的绘图系统，非常适合制作静态、交互式和动画的可视化。  \n# matplotlib.pyplot是matplotlib的绘图框架，提供了大量的绘图函数。  \nimport matplotlib.pyplot as plt  \nfrom matplotlib import animation, rc  # animation用于制作动画，rc用于配置matplotlib的全局设置  \n# 导入plotly.express库，Plotly Express是Plotly的Python图形库的高级封装，  \n# 提供了简单的API来创建复杂的交互式图表，非常适合数据可视化。  \nimport plotly.express as px  \n# 导入seaborn库，Seaborn是基于matplotlib的Python数据可视化库，  \n# 提供了一个高级接口来绘制吸引人的统计图形，如热图、箱线图、小提琴图等。  \nimport seaborn as sns  \n# 导入pydicom库，pydicom是一个用于处理DICOM文件的Python库，  \n# DICOM文件常用于医学影像的存储和传输，pydicom提供了读取、修改和写入DICOM文件的功能。  \nimport pydicom \n# 导入json库，json模块提供了一种很简单的方式来编码和解码JSON数据，  \n# JSON（JavaScript Object Notation）是一种轻量级的数据交换格式。  \nimport json  \n# 导入glob库，glob模块提供了一个函数用于从目录通配符搜索中生成文件列表，  \n# 这在处理具有特定模式的文件集合时非常有用。  \nimport glob  \n# 导入collections模块，collections模块包含了除列表、字典、集合、元组之外的特殊容器数据类型，  \n# 如Counter、defaultdict等，这些容器提供了额外的功能。  \nimport collections  \n# 从pydicom的pixel_data_handlers.util模块导入apply_voi_lut函数，  \n# 该函数用于应用灰度查找表（VOI LUT），这在处理医学影像的DICOM文件时非常有用，  \n# 可以将原始的像素值转换为更有意义的灰度值。  \nfrom pydicom.pixel_data_handlers.util import apply_voi_lut  \n# 导入cv2库，即OpenCV（Open Source Computer Vision Library），  \n# 是一个跨平台的计算机视觉库，提供了图像处理、视频分析、机器学习等功能。  \nimport cv2  \n# 导入skimage的measure模块，skimage（scikit-image）是Python的一个图像处理库，  \n# measure模块提供了图像处理中的测量功能，如连通区域检测、边缘检测等。  \nfrom skimage import measure  \n# 导入mpl_toolkits.mplot3d.art3d的Poly3DCollection，  \n# 这是matplotlib的一个工具，用于在3D图中绘制多边形，非常适合创建3D可视化。  \nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection  \n# 导入plotly.graph_objects库，Plotly Graph Objects是用于创建图表的类集合，  \n# 它提供了比Plotly Express更低层次的接口，允许更细粒度的控制图表的各个方面。  \nimport plotly.graph_objects as go  \n# 导入random库，random模块用于生成随机数，  \n# 在模拟、测试、游戏开发等领域非常有用。  \nimport random  \n# 导入mpl_toolkits.mplot3d的Axes3D，这是matplotlib中用于创建3D图形的类。  \nfrom mpl_toolkits.mplot3d import Axes3D  \n# 可选：过滤警告信息，忽略所有警告。  \n# import warnings  \n# warnings.filterwarnings('ignore')","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:40.914745Z","iopub.execute_input":"2024-07-24T04:09:40.915118Z","iopub.status.idle":"2024-07-24T04:09:40.925560Z","shell.execute_reply.started":"2024-07-24T04:09:40.915088Z","shell.execute_reply":"2024-07-24T04:09:40.924566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"二、读取数据","metadata":{}},{"cell_type":"code","source":"# 从指定路径读取训练集的标签坐标CSV文件到DataFrame label_coordinates_df中  \n# 这个文件包含了图像中特定标签（如腰椎退行性变化）的坐标信息  \nlabel_coordinates_df = pd.read_csv('/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_label_coordinates.csv')  \n# 从指定路径读取训练集序列描述CSV文件到DataFrame train_series中  \n# 这个文件包含了训练集中每个图像序列的描述信息  \ntrain_series = pd.read_csv('/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_series_descriptions.csv')  \n# 从指定路径读取训练集CSV文件到DataFrame df_train中  \n# 这个文件包含了训练集中每个图像的基本信息，如ID、标签等  \ndf_train = pd.read_csv('/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train.csv')  \n# 从指定路径读取样本提交CSV文件到DataFrame df_sub中  \n# 这个文件是一个模板，用于提交竞赛结果，包含了测试集中图像的ID和占位符的预测值  \ndf_sub = pd.read_csv('/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/sample_submission.csv')  \n# 从指定路径读取测试集序列描述CSV文件到DataFrame test_series中  \n# 这个文件包含了测试集中每个图像序列的描述信息，有助于理解测试数据  \ntest_series = pd.read_csv('/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/test_series_descriptions.csv')","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:40.927779Z","iopub.execute_input":"2024-07-24T04:09:40.928194Z","iopub.status.idle":"2024-07-24T04:09:41.035294Z","shell.execute_reply.started":"2024-07-24T04:09:40.928159Z","shell.execute_reply":"2024-07-24T04:09:41.034381Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"三、加载图像","metadata":{}},{"cell_type":"code","source":"# 设置一个文件夹路径，该路径指向训练图像中的一个特定研究（study）  \nfolder_path = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/100206310/1012284084'  \n# 从该文件夹中列出所有以'.dcm'结尾的文件，即DICOM图像文件  \ndicom_files = [f for f in os.listdir(folder_path) if f.endswith('.dcm')]  \n# 读取包含标签坐标的CSV文件  \nlabel_coordinates_df = pd.read_csv('/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_label_coordinates.csv')  \n# 从文件夹路径中提取研究ID（study_id），假设路径中的倒数第二个元素是研究ID  \nstudy_id = folder_path.split('/')[-2]  \n# 从标签坐标DataFrame中筛选出与当前研究ID相对应的行  \nstudy_label_coordinates = label_coordinates_df[label_coordinates_df['study_id'] == int(study_id)]  \n# 初始化两个空列表，用于存储有对应标签坐标的DICOM文件名和相应的标签坐标  \nfiltered_dicom_files = []  \nfiltered_label_coordinates = []  \n# 遍历所有DICOM文件  \nfor dicom_file in dicom_files:  \n    # 从文件名中提取实例编号（instance_number），假设文件名（不含扩展名）即为实例编号  \n    instance_number = int(dicom_file.split('.')[0])  \n    # 在标签坐标DataFrame中查找与当前实例编号相匹配的行  \n    corresponding_coordinates = study_label_coordinates[study_label_coordinates['instance_number'] == instance_number]  \n    # 如果找到了匹配的标签坐标，则将该DICOM文件名和标签坐标分别添加到两个列表中  \n    if not corresponding_coordinates.empty:  \n        filtered_dicom_files.append(dicom_file)  \n        filtered_label_coordinates.append(corresponding_coordinates)  \n# 使用matplotlib绘制图像  \n# 这里假设我们想要绘制第二行（索引从1开始）的四个DICOM图像及其标签坐标  \nfig, axs = plt.subplots(1, 4, figsize=(20, 5))  \nsecond_row_index = 1\nsecond_row_images = filtered_dicom_files[second_row_index:second_row_index+4]  \nsecond_row_coordinates = filtered_label_coordinates[second_row_index:second_row_index+4]  \n# 遍历选定的DICOM图像和对应的标签坐标  \nfor i, (dicom_file, label_coordinates) in enumerate(zip(second_row_images, second_row_coordinates)):  \n    # 构造DICOM文件的完整路径  \n    dicom_file_path = os.path.join(folder_path, dicom_file)  \n    # 使用pydicom库读取DICOM文件  \n    dicom_data = pydicom.dcmread(dicom_file_path)  \n    # 获取图像数据  \n    image = dicom_data.pixel_array  \n    # 在matplotlib的子图上显示图像  \n    axs[i].imshow(image, cmap='gray')  \n    axs[i].set_title(f'DICOM Image - {dicom_file}')  \n    axs[i].axis('off')  # 关闭坐标轴  \n    # 遍历标签坐标，并在图像上绘制标记  \n    for _, row in label_coordinates.iterrows():  \n        axs[i].plot(row['x'], row['y'], 'ro', markersize=5)  # 假设'x'和'y'是标签坐标的列名  \n# 调整子图布局以避免重叠  \nplt.tight_layout()    \n# 显示图像  \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:41.062560Z","iopub.execute_input":"2024-07-24T04:09:41.062831Z","iopub.status.idle":"2024-07-24T04:09:42.410383Z","shell.execute_reply.started":"2024-07-24T04:09:41.062809Z","shell.execute_reply":"2024-07-24T04:09:42.409482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_dicom(path):  \n    \"\"\"  \n    加载DICOM文件，对其像素数据进行预处理，并返回处理后的数据。  \n    参数:  \n    path (str): DICOM文件的路径。  \n    返回:  \n    numpy.ndarray: 预处理后的像素数据，数据类型为np.uint8，且值域为0-255。  \n    过程说明:  \n    1. 使用pydicom库的read_file函数读取指定路径的DICOM文件。  \n    2. 从DICOM文件中获取像素数组（pixel_array），即图像的原始数据。  \n    3. 对像素数组进行初步处理，通过减去数组中的最小值，确保所有像素值均为非负，这一步是简单的归一化预处理。  \n    4. 检查处理后的像素数组的最大值是否不为0，这是为了避免在后续步骤中进行除以零的操作。  \n    5. 如果最大值不为0，则将像素数组的每个元素除以最大值，完成归一化处理，使得像素值的范围变为0到1之间。  \n    6. 将归一化后的像素值乘以255，这样做的目的是将像素值的范围从0到1转换为0到255，这是大多数图像处理库和显示设备所期望的数值范围。  \n    7. 使用astype方法将处理后的像素数组的数据类型转换为np.uint8，这是图像数据常用的数据类型，因为它可以高效地存储和传输图像信息。  \n    8. 返回处理后的像素数组，供进一步分析或显示使用。  \n    \"\"\"  \n    dicom = pydicom.read_file(path)  # 读取DICOM文件  \n    data = dicom.pixel_array  # 获取DICOM文件中的像素数组  \n    data = data - np.min(data)  # 减去数组中的最小值，确保没有负值  \n    if np.max(data) != 0:  # 检查最大值是否不为0，以避免除以零  \n        data = data / np.max(data)  # 除以最大值进行归一化  \n    data = (data * 255).astype(np.uint8)  # 缩放值域并转换为np.uint8类型  \n    return data  # 返回处理后的像素数组","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:42.440416Z","iopub.execute_input":"2024-07-24T04:09:42.440710Z","iopub.status.idle":"2024-07-24T04:09:42.447911Z","shell.execute_reply.started":"2024-07-24T04:09:42.440686Z","shell.execute_reply":"2024-07-24T04:09:42.447118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"四、图像动画视图","metadata":{}},{"cell_type":"code","source":"# 设置matplotlib的动画输出格式为HTML，这样可以在Jupyter Notebook中直接显示动画  \nrc('animation', html='jshtml')  \n# 定义一个函数来加载单个DICOM文件  \ndef load_dicom(filename):  \n    try:  \n        # 使用pydicom库读取DICOM文件  \n        ds = pydicom.dcmread(filename)  \n        # 返回图像数据  \n        return ds.pixel_array  \n    except Exception as e:  \n        # 如果读取文件时发生错误，打印错误信息并返回None  \n        print(f\"Error reading DICOM file {filename}: {e}\")  \n        return None  \n# 定义一个函数来加载文件夹内的一系列DICOM文件  \ndef load_dicom_line(path):  \n    # 使用glob模块找到路径下的所有文件，并按文件名中的数字排序  \n    t_paths = sorted(  \n        glob.glob(os.path.join(path, \"*\")),  \n        key=lambda x: int(os.path.splitext(os.path.basename(x))[0].split(\"-\")[-1]),  \n    )  \n    images = []  \n    for filename in t_paths:  \n        # 加载每个DICOM文件  \n        data = load_dicom(filename)  \n        # 如果数据有效且不是全黑（最大值大于0），则添加到列表中  \n        if data is not None and data.max() > 0:  \n            images.append(data)  \n    return images  ","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:42.449122Z","iopub.execute_input":"2024-07-24T04:09:42.449478Z","iopub.status.idle":"2024-07-24T04:09:42.803462Z","shell.execute_reply.started":"2024-07-24T04:09:42.449453Z","shell.execute_reply":"2024-07-24T04:09:42.802538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 定义一个函数来创建动画  \ndef create_animation(ims):  \n    # 创建一个图形窗口  \n    fig = plt.figure(figsize=(6, 6))  \n    # 关闭坐标轴显示  \n    plt.axis('off')  \n    # 显示第一张图像  \n    im = plt.imshow(ims[0], cmap=\"gray\")  \n    # 定义一个更新图像的函数  \n    def animate_func(i):  \n        im.set_array(ims[i])  # 更新图像数据  \n        return [im]  # 返回需要更新的艺术家列表  \n    # 创建动画  \n    return animation.FuncAnimation(fig, animate_func, frames=len(ims), interval=1000//24)  \n# 指定包含DICOM文件的文件夹路径  \npath_to_folder = \"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/100206310/1792451510\"  \n# 加载文件夹内的DICOM图像  \nimages = load_dicom_line(path_to_folder)  \n# 如果找到了有效的图像，则创建动画并显示  \nif images:  \n    anim = create_animation(images)  \n    # 在Jupyter Notebook中显示动画  \n    plt.show()  \nelse:  \n    # 如果没有找到有效的图像，则打印消息  \n    print(\"No valid images found.\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"五、绘制直方图","metadata":{}},{"cell_type":"code","source":"def read_dicom_files(path_to_folder, num_files=5):  \n    # 构造一个文件模式，用于匹配指定文件夹内所有的.dcm（DICOM）文件  \n    files_glob = os.path.join(path_to_folder, \"*.dcm\")  \n    # 使用glob.glob函数查找所有匹配的DICOM文件，并通过排序确保文件按特定顺序（基于文件名中的数字）返回  \n    # 这里假设文件名中包含可以表示顺序的数字，并且这些数字是文件名中'-'分隔的最后一个部分  \n    dicom_files = sorted(glob.glob(files_glob), key=lambda f: int(os.path.splitext(os.path.basename(f))[0].split('-')[-1]))  \n    # 使用pydicom.dcmread函数读取排序后的前num_files个DICOM文件，并返回这些文件的DICOMDataset对象列表  \n    return [pydicom.dcmread(f) for f in dicom_files[:num_files]]","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:42.804870Z","iopub.execute_input":"2024-07-24T04:09:42.805494Z","iopub.status.idle":"2024-07-24T04:09:42.812718Z","shell.execute_reply.started":"2024-07-24T04:09:42.805459Z","shell.execute_reply":"2024-07-24T04:09:42.811664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def calculate_level(mean, std):  \n    # 根据给定的均值（mean）和标准差（std），计算一个阈值（level）  \n    # 这个阈值是基于均值加上1.7倍的标准差计算得出的，通常用于图像处理中的阈值设置  \n    return mean + 1.7 * std","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:42.817413Z","iopub.execute_input":"2024-07-24T04:09:42.817708Z","iopub.status.idle":"2024-07-24T04:09:42.826670Z","shell.execute_reply.started":"2024-07-24T04:09:42.817673Z","shell.execute_reply":"2024-07-24T04:09:42.825721Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def stats_image(image):  \n    # 计算图像中非零像素的均值和标准差  \n    # 首先，找到所有非零像素  \n    noncero_pixels = image[np.nonzero(image)]  \n    # 如果图像中没有非零像素（即全黑图像），则均值和标准差都设为0  \n    if noncero_pixels.size == 0:  \n        mean = 0  \n        std = 0  \n    else:  \n        # 否则，计算非零像素的均值和标准差  \n        mean = np.mean(noncero_pixels)  \n        std = np.std(noncero_pixels)  \n    # 返回均值和标准差  \n    return mean, std","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:42.827855Z","iopub.execute_input":"2024-07-24T04:09:42.828181Z","iopub.status.idle":"2024-07-24T04:09:42.836768Z","shell.execute_reply.started":"2024-07-24T04:09:42.828149Z","shell.execute_reply":"2024-07-24T04:09:42.835844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def image_orientation(dicom):  \n    # 根据DICOM文件中的ImageOrientationPatient字段判断图像的朝向  \n    # 首先，假设朝向未知  \n    rt = 'unknown'  \n    # 提取ImageOrientationPatient字段的值，并四舍五入到整数 \n    x1, y1, _, x2, y2, _ = [round(v) for v in dicom.ImageOrientationPatient]  \n    # 根据ImageOrientationPatient的值判断图像的朝向  \n    # 这里只考虑了三种常见的朝向：冠状面（coronal）、轴向面（axial）和矢状面（sagittal）  \n    if (x1, y1, x2, y2) == (1, 0, 0, 0):  \n        rt = 'coronal'  \n    elif (x1, y1, x2, y2) == (1, 0, 0, 1):  \n        rt = 'axial'  \n    elif (x1, y1, x2, y2) == (0, 1, 0, 0):  \n        rt = 'sagittal'  \n    # 返回图像的朝向  \n    return rt","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:42.838254Z","iopub.execute_input":"2024-07-24T04:09:42.838597Z","iopub.status.idle":"2024-07-24T04:09:42.846598Z","shell.execute_reply.started":"2024-07-24T04:09:42.838566Z","shell.execute_reply":"2024-07-24T04:09:42.845747Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 定义一个函数来绘制图像的直方图  \ndef plot_image_hist(image):  \n    # 计算图像的均值和标准差  \n    mean, std = stats_image(image)    \n    # 将图像数据展平为一维数组  \n    pixels = image.ravel()  \n    # 过滤出非零像素  \n    noncero_pixels = pixels[np.nonzero(pixels)]\n    # 将非零像素标准化（去均值并除以标准差）  \n    noncero_pixels = (noncero_pixels - mean) / std  \n    # 计算并获取超过某个阈值（基于均值和标准差计算）的非零像素数量  \n    over_threshold = np.count_nonzero(noncero_pixels > calculate_level(mean, std))  \n    # 创建包含两个子图的图形  \n    fig, (axi, axh) = plt.subplots(1, 2, figsize=(20, 3), gridspec_kw={'width_ratios': [1, 4]})  \n    # 设置图形的总标题，包括超过阈值的非零像素数量  \n    fig.suptitle(f'scan # ({over_threshold})')  \n    # 在右侧子图上绘制非零像素的直方图  \n    axh.hist(noncero_pixels, 200, range=(-5, 5))  # 绘制直方图，范围限制在-5到5之间  \n    axh.set_xlim(-5, 5)  # 设置x轴限制  \n    # 获取直方图的y轴限制  \n    ax_limits = axh.get_ylim()  \n    # 在直方图上绘制均值、均值+标准差和阈值的垂直线  \n    axh.vlines(mean, ymin=ax_limits[0], ymax=ax_limits[1], colors='r', label='Mean')  \n    axh.vlines(mean + std, ymin=ax_limits[0], ymax=ax_limits[1], colors='g', linestyles='dotted', label='Mean + Std')  \n    axh.vlines(calculate_level(mean, std), ymin=ax_limits[0], ymax=ax_limits[1], colors='b', linestyles='dashed', label='Threshold')  \n    # 在左侧子图上显示原始图像  \n    axi.imshow(image, cmap=plt.cm.gray)  # 使用灰度图显示  \n    axi.grid(False)  # 不显示网格  \n    axi.axis('off')  # 关闭坐标轴  \n    axh.legend()  # 显示直方图上的图例  \n    plt.show()  # 显示图形  \n# 指定DICOM文件的路径  \npath_to_folder = \"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/1002894806/866293114\"  \n# read_dicom_files是一个自定义函数，用于读取指定文件夹中的所有DICOM文件  \ndicom_files = read_dicom_files(path_to_folder)  \n# 如果找到DICOM文件，则对第一张图像进行可视化  \nif dicom_files:  \n    first_image = dicom_files[0].pixel_array  # 获取第一张图像的像素数组  \n    plot_image_hist(first_image)  # 调用函数进行可视化  \nelse:  \n    print(\"No DICOM files found.\")  # 如果没有找到DICOM文件，则打印消息","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:42.847829Z","iopub.execute_input":"2024-07-24T04:09:42.848195Z","iopub.status.idle":"2024-07-24T04:09:43.587397Z","shell.execute_reply.started":"2024-07-24T04:09:42.848160Z","shell.execute_reply":"2024-07-24T04:09:43.586456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 定义一个函数来读取指定文件夹中的DICOM文件  \n# num_files参数控制读取的文件数量  \ndef read_dicom_files(path_to_folder, num_files=5):  \n    # 使用glob模式匹配文件夹中的所有.dcm文件  \n    files_glob = os.path.join(path_to_folder, \"*.dcm\")  \n    # 对匹配到的文件进行排序，这里假设文件名中包含可以按数字排序的序列信息  \n    dicom_files = sorted(glob.glob(files_glob), key=lambda f: int(os.path.splitext(os.path.basename(f))[0].split('-')[-1]))  \n    # 读取前num_files个DICOM文件  \n    return [pydicom.dcmread(f) for f in dicom_files[:num_files]]  \n# 定义一个函数来获取DICOM文件中的像素数组  \ndef get_flair_images(dicom_files):  \n    # 遍历DICOM文件，提取像素数组  \n    images = [s.pixel_array for s in dicom_files]  \n    # 将像素数组列表转换为NumPy数组  \n    return np.array(images)  \n# 指定DICOM文件的路径  \npath_to_folder = \"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/100206310/1012284084\"  \n# 读取DICOM文件  \ndicom_files = read_dicom_files(path_to_folder)  \n# 获取FLAIR图像的像素数组  \nflair_images = get_flair_images(dicom_files)  \nif flair_images.ndim == 3:  # 确保是三维数组（多切片）  \n    first_frame = flair_images[0, :, :]  # 提取第一帧  \n    fig = px.imshow(first_frame, labels=dict(x=\"Pixel Column\", y=\"Pixel Row\"))  \n    fig.update_layout(title=\"First Frame of FLAIR Image\")  \n    fig.show()  ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-07-24T04:09:43.588680Z","iopub.execute_input":"2024-07-24T04:09:43.588980Z","iopub.status.idle":"2024-07-24T04:09:43.701237Z","shell.execute_reply.started":"2024-07-24T04:09:43.588955Z","shell.execute_reply":"2024-07-24T04:09:43.700295Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 遍历flair_images数组中的前5个图像（如果存在的话）  \nfor img in flair_images[:5]:  \n    # 对每个图像调用plot_image_hist函数，可能用于绘制图像的直方图   \n    # 并以某种方式（可能是绘图）显示该图像的直方图  \n    plot_image_hist(img)","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:43.702640Z","iopub.execute_input":"2024-07-24T04:09:43.702949Z","iopub.status.idle":"2024-07-24T04:09:47.059238Z","shell.execute_reply.started":"2024-07-24T04:09:43.702921Z","shell.execute_reply":"2024-07-24T04:09:47.058420Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"六、图像的信息","metadata":{}},{"cell_type":"code","source":"# 获取第二个DICOM文件的路径  \ndicom_files[1]","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:47.060331Z","iopub.execute_input":"2024-07-24T04:09:47.060600Z","iopub.status.idle":"2024-07-24T04:09:47.068561Z","shell.execute_reply.started":"2024-07-24T04:09:47.060575Z","shell.execute_reply":"2024-07-24T04:09:47.067591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"七、图像的3D视图","metadata":{}},{"cell_type":"code","source":"# 定义函数来加载DICOM图像  \ndef load_dicom_images(folder_path, num_images=5):  \n    # 遍历指定文件夹下的所有文件，筛选出以.dcm结尾的DICOM文件  \n    dicom_files = sorted([os.path.join(folder_path, f) for f in os.listdir(folder_path) if f.endswith('.dcm')])  \n    # 只取前num_images个文件（如果文件夹中文件多于num_images）  \n    dicom_files = dicom_files[:num_images]  \n    # 初始化一个空列表来存储图像数据  \n    images = []  \n    # 遍历DICOM文件列表，读取每个文件的像素数组并添加到列表中  \n    for f in dicom_files:  \n        img = pydicom.dcmread(f).pixel_array  \n        images.append(img)  \n    # 使用numpy的stack函数将图像列表沿第一个轴（axis=0）堆叠成一个3D数组  \n    # 这里假设每个图像切片在空间上是连续的，并且我们想要沿Z轴堆叠它们  \n    # 因此，堆叠后的数组形状将是(num_images, height, width)或(num_images, width, height)，取决于图像的原始尺寸  \n    return np.stack(images, axis=0)  \n# 定义函数来绘制3D图像中的等值面  \ndef plot_3d(image, threshold=-300):  \n    # 使用marching_cubes算法从3D图像中提取等值面  \n    # 这里，threshold是确定等值面位置的阈值  \n    verts, faces, _, _ = marching_cubes(image, level=threshold)  \n    # 创建一个新的图形和3D子图  \n    fig = plt.figure(figsize=(10, 10))  \n    ax = fig.add_subplot(111, projection='3d')  \n    # 创建一个Poly3DCollection对象来绘制多边形集合  \n    # verts[faces]给出了等值面的顶点坐标  \n    mesh = Poly3DCollection(verts[faces], alpha=0.1)  \n    # 设置多边形的面颜色  \n    face_color = [0.5, 0.5, 1]  # 浅蓝色  \n    mesh.set_facecolor(face_color)  \n    # 将多边形集合添加到3D子图中  \n    ax.add_collection3d(mesh)  \n    # 设置3D子图的x、y、z轴的限制，以匹配图像的尺寸  \n    ax.set_xlim(0, image.shape[0])  \n    ax.set_ylim(0, image.shape[1])  \n    ax.set_zlim(0, image.shape[2])  \n    # 显示图形  \n    plt.show()  \n# 指定DICOM图像文件夹的路径  \nfolder_path = \"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/100206310/1012284084\"  \n# 加载DICOM图像并堆叠成3D数组  \ndicom_images = load_dicom_images(folder_path, num_images=5)  \n# 绘制3D图像中的等值面，这里使用了一个示例阈值300  \n# 注意：这个阈值可能需要根据您的具体数据集进行调整  \nplot_3d(dicom_images, threshold=300)","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:09:47.070527Z","iopub.execute_input":"2024-07-24T04:09:47.070992Z","iopub.status.idle":"2024-07-24T04:10:03.367088Z","shell.execute_reply.started":"2024-07-24T04:09:47.070954Z","shell.execute_reply":"2024-07-24T04:10:03.366164Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 定义函数来加载DICOM图像  \ndef load_dicom_images(folder_path, num_images=5):  \n#     加载指定文件夹下的DICOM图像文件，并返回它们作为numpy数组的堆叠。  \n#     参数:  \n#     folder_path (str): 包含DICOM文件的文件夹路径。  \n#     num_images (int): 要加载的DICOM图像数量，默认为5。  \n#     返回:  \n#     numpy.ndarray: 一个3D numpy数组，其中包含了堆叠的DICOM图像切片。  \n#                    如果图像尺寸不一致，这可能会导致问题。  \n    # 列出文件夹中所有以.dcm结尾的文件，并排序  \n    # 然后取前num_images个文件  \n    dicom_files = sorted([os.path.join(folder_path, f) for f in os.listdir(folder_path) if f.endswith('.dcm')])[:num_images]  \n    # 读取每个DICOM文件的像素数组，并将它们存储在一个列表中  \n    images = [pydicom.dcmread(f).pixel_array for f in dicom_files]  \n    # 注意：这里使用axis=-1来沿最后一个轴（即深度/切片方向）堆叠图像  \n    # 这假设所有图像切片都有相同的宽度和高度  \n    return np.stack(images, axis=-1)  \n# 定义函数来绘制交互式3D DICOM图像  \ndef plot_3d_interactive(image, threshold=-300):    \n#     使用Plotly库绘制3D DICOM图像的交互式可视化。  \n#     参数:  \n#     image (numpy.ndarray): 3D numpy数组，包含DICOM图像切片。  \n#     threshold (int): 用于marching_cubes算法的阈值，用于确定等值面的位置。  \n#     返回:  \n#     None: 直接显示交互式图表。  \n    # 使用marching_cubes算法从3D图像中提取等值面  \n    verts, faces, _, _ = measure.marching_cubes(image, level=threshold)  \n    # 创建一个Plotly Mesh3d对象，用于绘制3D网格  \n    fig = go.Figure(data=[  \n        go.Mesh3d(  \n            x=verts[:, 0],  # x坐标  \n            y=verts[:, 1],  # y坐标  \n            z=verts[:, 2],  # z坐标  \n            i=faces[:, 0],  # 网格顶点索引（三角形顶点的第一个顶点）  \n            j=faces[:, 1],  # 网格顶点索引（三角形顶点的第二个顶点）  \n            k=faces[:, 2],  # 网格顶点索引（三角形顶点的第三个顶点）  \n            color='blue',   # 网格颜色  \n            opacity=0.1     # 网格透明度  \n        )  \n    ])  \n    # 更新图表的布局设置  \n    fig.update_layout(  \n        scene=dict(  \n            xaxis=dict(visible=True),  # x轴可见  \n            yaxis=dict(visible=True),  # y轴可见  \n            zaxis=dict(visible=True)   # z轴可见  \n        ),  \n        width=800,  # 图表宽度  \n        height=800, # 图表高度  \n        title=\"Interactive 3D DICOM Image Visualization\"  # 图表标题  \n    )  \n    # 显示图表  \n    fig.show()  \n# 指定DICOM图像文件夹的路径  \nfolder_path = \"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/1002894806/866293114\"  \n# 加载DICOM图像并堆叠成3D数组  \ndicom_images = load_dicom_images(folder_path, num_images=5)  \n# 注意：这里使用的threshold值（100）可能需要根据实际图像进行调整  \n# 以确保能够正确提取出有意义的等值面  \nplot_3d_interactive(dicom_images, threshold=100)","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:10:03.368306Z","iopub.execute_input":"2024-07-24T04:10:03.368592Z","iopub.status.idle":"2024-07-24T04:10:03.613682Z","shell.execute_reply.started":"2024-07-24T04:10:03.368567Z","shell.execute_reply":"2024-07-24T04:10:03.612672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_3d_image(image):  \n#     绘制图像的三维表示。  \n#     参数:  \n#     image (numpy.ndarray): 要绘制的二维图像数组。  \n#     该函数使用matplotlib的3D绘图功能来生成图像的三维曲面图，其中X和Y轴代表图像的二维空间，  \n#     Z轴代表图像的强度（像素值）。  \n    fig = plt.figure(figsize=(10, 8))  # 创建一个大小为10x8英寸的图形窗口  \n    ax = fig.add_subplot(111, projection='3d')  # 添加一个3D子图  \n    rows, cols = image.shape  # 获取图像的行数和列数  \n    x, y = np.meshgrid(np.arange(cols), np.arange(rows))  # 创建网格坐标矩阵  \n    # 使用plot_surface函数绘制图像的三维曲面，cmap设置颜色映射，edgecolor设置为'none'以不显示边缘线  \n    ax.plot_surface(x, y, image, cmap='viridis', edgecolor='none')  \n    # 设置3D图的X、Y、Z轴标签和标题  \n    ax.set_xlabel('X')  \n    ax.set_ylabel('Y')  \n    ax.set_zlabel('Intensity')  \n    ax.set_title('3D Plot of Image')  \n    plt.show()  # 显示图形  \n# 假设flair_images是一个包含多个二维图像数组的列表或数组  \n# 这里我们绘制列表中的第一个图像  \nplot_3d_image(flair_images[0])  # 调用函数，传入flair_images列表中的第一个图像进行绘制","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:10:03.614910Z","iopub.execute_input":"2024-07-24T04:10:03.615194Z","iopub.status.idle":"2024-07-24T04:10:04.282172Z","shell.execute_reply.started":"2024-07-24T04:10:03.615168Z","shell.execute_reply":"2024-07-24T04:10:04.281273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"简单的数据准备训练/预测","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import LabelEncoder  # 用于将标签转换为整数编码  \nfrom sklearn.ensemble import RandomForestClassifier  # 随机森林分类器，虽然在这段代码中未直接使用  \n# df_train是已经存在的DataFrame，包含研究ID和多个条件级别的严重程度数据  \n# 使用melt函数将df_train中的非id_vars列（这里是除了study_id以外的所有列）转换为长格式数据  \n# var_name指定了变量名的列名，value_name指定了值的列名  \ndf_train_melted = df_train.melt(id_vars=['study_id'], var_name='condition_level', value_name='severity')  \n# 将condition_level列拆分为'condition'和'level'两列  \n# 使用str.rsplit('_', n=1, expand=True)方法，以'_'为分隔符，从右边开始拆分，并扩展为两列  \ndf_train_melted[['condition', 'level']] = df_train_melted['condition_level'].str.rsplit('_', n=1, expand=True)  \n# 使用LabelEncoder对'severity'列进行编码，将标签转换为整数  \nle_severity = LabelEncoder()  \ndf_train_melted['severity_encoded'] = le_severity.fit_transform(df_train_melted['severity'])  \n# 准备训练数据  \n# 选取'study_id', 'condition', 'level'作为特征  \nX_train = df_train_melted[['study_id', 'condition', 'level']]  \n# 'severity_encoded'作为目标变量  \ny_train = df_train_melted['severity_encoded']  \n# 对特征进行独热编码（One-Hot Encoding），以处理分类变量  \nX_train = pd.get_dummies(X_train, columns=['condition', 'level'])  \n# 准备测试数据  \n# 假设test_series是包含测试数据的Series或DataFrame，这里我们假设它至少包含'study_id'  \ntest_rows = []  \n# 遍历test_series中的每一行  \nfor _, row in test_series.iterrows():  \n    # 遍历所有可能的条件和级别组合  \n    for condition in ['left_neural_foraminal_narrowing', 'right_neural_foraminal_narrowing', 'left_subarticular_stenosis', 'right_subarticular_stenosis', 'spinal_canal_stenosis']:  \n        for level in ['l1_l2', 'l2_l3', 'l3_l4', 'l4_l5', 'l5_s1']:  \n            # 为每个条件和级别的组合创建一个测试行  \n            test_rows.append({  \n                'study_id': row['study_id'],  \n                'condition': condition,  \n                'level': level  \n            })  \n# 将test_rows转换为DataFrame  \nX_test = pd.DataFrame(test_rows)  \n# 对X_test进行独热编码  \nX_test = pd.get_dummies(X_test, columns=['condition', 'level'])  \n# 确保X_test的列与X_train的列一致，如果X_test中有X_train中没有的列，则这些列将被丢弃  \n# 如果X_train中有X_test中没有的列，则这些列将在X_test中以0填充  \nX_test = X_test.reindex(columns=X_train.columns, fill_value=0)","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:10:04.339960Z","iopub.execute_input":"2024-07-24T04:10:04.340217Z","iopub.status.idle":"2024-07-24T04:10:04.471822Z","shell.execute_reply.started":"2024-07-24T04:10:04.340194Z","shell.execute_reply":"2024-07-24T04:10:04.470971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 初始化一个RandomForestClassifier模型实例  \n# RandomForestClassifier是随机森林分类器，它基于决策树算法构建了一个由多个决策树组成的森林  \n# 这些决策树在训练数据集的不同子集上独立训练，并通过平均它们的预测结果来改进整体预测的准确性  \n# 这里没有指定任何参数，因此将使用RandomForestClassifier的默认参数设置  \nmodel = RandomForestClassifier()  \n# 使用训练数据X_train和训练标签y_train来拟合（训练）模型  \n# 这个过程涉及到在X_train提供的特征上训练随机森林中的每一棵树  \n# 每棵树都会尝试找到最佳的特征和分割点，以便将数据集分割成越来越纯净的子集  \n# 最终，整个随机森林会基于所有树的预测结果来做出最终的分类决策  \nmodel.fit(X_train, y_train)  \n# 经过fit方法后，model对象现在包含了训练好的随机森林模型  ","metadata":{"execution":{"iopub.status.busy":"2024-07-24T04:10:04.504786Z","iopub.execute_input":"2024-07-24T04:10:04.505032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 使用训练好的随机森林模型对测试集X_test进行预测，但这次不是直接预测类别标签，而是预测每个类别的概率  \n# predict_proba方法返回的是一个二维数组，其中每一行对应X_test中的一个样本，每一列对应一个类别的预测概率  \npredictions_proba = model.predict_proba(X_test)  \n# 将预测的概率转换为DataFrame格式，以便更方便地处理和查看结果  \n# DataFrame的列名被设置为le_severity.classes_，这通常是使用LabelEncoder或类似工具对目标变量y进行编码时得到的类别名称  \n# 这样做可以确保DataFrame的列与预测概率的顺序相匹配  \npredictions_df = pd.DataFrame(predictions_proba, columns=le_severity.classes_)  \n# 将测试集X_test中的study_id列添加到predictions_df中，以便将预测结果与原始样本关联起来  \npredictions_df['study_id'] = X_test['study_id'].values  \n# 创建一个新的列condition_level，它通过将测试集中的索引映射到test_rows DataFrame的相应行，并格式化condition和level列来得到  \n# 这里假设test_rows是一个包含测试集原始数据的DataFrame，且包含condition和level列  \n# index.map函数用于将X_test的索引（即样本在原始测试集中的位置）映射到test_rows DataFrame的对应行上，并提取所需的信息  \n# 最后，通过字符串格式化将condition和level组合成一个新的字符串  \npredictions_df['condition_level'] = X_test.index.map(lambda idx: f\"{test_rows[idx]['condition']}_{test_rows[idx]['level']}\")  \n# 创建一个新的列row_id，它是通过将study_id和condition_level组合成一个唯一的字符串来得到的  \n# 这个字符串用于在提交预测结果时作为行的唯一标识符  \n# 注意：这里study_id是字符串类型，如果不是，则通过astype(str)将其转换为字符串类型  \npredictions_df['row_id'] = predictions_df['study_id'].astype(str) + '_' + predictions_df['condition_level']  \n# 这段代码的目的是为了准备预测结果的DataFrame，以便后续可以将其保存为CSV文件或用于其他分析  \n# 预测结果DataFrame包含了每个测试样本的预测概率、study_id、condition_level以及用于标识每行的row_id","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 使用pandas库的read_csv函数来读取一个CSV文件\ndf_sub1 = pd.read_csv('/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/sample_submission.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 从predictions_df DataFrame中获取'Normal/Mild'列的第一行值，并将其存储在normal_mild_value变量中  \nnormal_mild_value = (predictions_df['Normal/Mild'].iloc[0])  \n# 从predictions_df DataFrame中获取'Moderate'列的第一行值，并将其存储在moderate_value变量中  \nmoderate_value = (predictions_df['Moderate'].iloc[0])  \n# 从predictions_df DataFrame中获取'Severe'列的第一行值，并将其存储在severe_value变量中  \nsevere_value = (predictions_df['Severe'].iloc[0])  \n# 将'Normal/Mild'的预测值除以2.3，然后将结果存储到df_sub DataFrame的'normal_mild'列中  \n# 这里可能是为了对预测值进行某种形式的标准化或调整  \ndf_sub['normal_mild'] = normal_mild_value / 2.3  \n# 将'Moderate'的预测值乘以1.47，然后将结果存储到df_sub DataFrame的'moderate'列中  \ndf_sub['moderate'] = moderate_value * 1.47  \n# 计算'severe'列的值  \ndf_sub['severe'] = (1 - (normal_mild_value / 2.3 + moderate_value * 1.47)) + severe_value - severe_value  ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 将df_sub DataFrame保存到CSV文件中，文件名为'submission.csv'，不包括索引列  \ndf_sub.to_csv('submission.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"参考文献：\nhttps://www.kaggle.com/code/satyaprakashshukl/rsna-lumbar-spine-analysis","metadata":{}}]}