【2】Kaggle 医学影像数据读取

赛题称号：RSNA 2024 Lumbar Spine Degenerative Classification
中文：腰椎退行性病变分类

kaggle官网赛题链接：https://www.kaggle.com/competitions/rsna-2024-lumbar-spine-degenerative-classification/overview

文章安排

①、如何用python读取dcm/dicom文件
②、根据matplotlib可视化
③、制造频率散布直方图
④、代码汇总

文件依靠

# requirements.txt
# Python version 3.11.8
torch==2.3.1
torchvision==0.18.1
matplotlib==3.8.4
pydicom==2.4.4
numpy==1.26.4

pip install -r requirements.txt

读取dicom图画并做预处理

概述

本文中采纳pydicom包读取dicom文件，其要害代码格局为：

dcm_tensor = pydicom.dcmread(dcm_file)

留意数据集的途径，其在train_images文件下寄存了每一患者的数据，关于每一患者包含三张MRI图画，每张MRI图画寄存为一个文件夹。
需求留意的是，MRI图画为三维图画（dicom格局），一般习惯性将其每个切片别离保存为一个dcm文件，因而一张dicom图画将被存为一个文件夹，如下图

咱们能够选用如下途径拜访该dicom文件：

"./train_images/4003253/702807833"

读取途径

为了读取dicom图画，咱们需求写代码读取文件夹中的一切dcm文件

# dicom文件途径
dicom_dir = "./train_images/4003253/702807833"
# 保存一切dcm文件的途径
dicom_files = [os.path.join(dicom_dir, f) for f in os.listdir(dicom_dir) if f.endswith('.dcm')] 

• os.listdir：回来dicom_dir途径下的一切文件
• f.endswith('.dcm') ：挑选一切dcm格局的文件
• os.path.join: 将dcm文件名添加到dicom_dir之后
  暗示："./hello"+“1.dcm”->"./hello/1.dcm"

途径排序

这次的kaggle赛题所给的数据会集，文件名的迭代办法为：

1.dcm、2.dcm、...、9.dcm、10.dcm、11.dcm、...

这给咱们带来了必定的费事，由于在os的文件名排序规矩中，首要检索高位字母的ASCII码巨细做排序，也就是说10.dcm将被认为是2.dcm前面的文件。
对此，本文选用正则表达式的办法，完成了根据文件名中数字巨细排序。

def extract_number(filepath):
    # 获取文件名（包含扩展名）
    filename = os.path.basename(filepath)
    # 提取文件名中的数字部分，假定文件名以数字结束，如 '1.dcm'
    match = re.search(r'(\d+)\.dcm$', filename)
    return int(match.group(1)) if match else float('inf')

# 根据数字句柄排序
dicom_files.sort(key=extract_number)

该代码作用如下：

读取图画

为读取dicom图画，咱们需求顺次读取每一个dcm文件，并将其终究打包为3D tensor，下述代码完成了该功用：

# 创立空列表保存一切dcm文件
dcm_list= []

# 迭代每一个文件
for dcm_file in dicom_files:
    # 读取文件
    dcm = pydicom.dcmread(dcm_file)
    # 将其转为numpy格局
    image_data = dcm.pixel_array.astype(np.float32)
    # 参加文件列表 
    dcm_list.append(image_data)

# 将图片堆叠为3D张量
tensor_dcm = torch.stack([torch.tensor(image_data) for image_data in dcm_list])

数据预处理

常见的预处理办法有两种，归一化（Normalization）和量化（Quantization）

• 归一化：将数据缩放到某个规范规模内的进程。常见的归一化办法包含最小-最大归一化（Min-Max Normalization）和Z-score规范化（Z-score Normalization），前者将数据归一化至[0,1]规模，后者将数据转化为规范正态散布。本例中选用Min-Max计划。

• 量化：量化是将数据的值域退化到离散值的进程。常用于削减存储和计算成本，尤其在神经网络模型中。量化通常将浮点数值转换为整数值。量化前一般先进行归一化。

归一化的完成如下：

def norm_tensor(tensor_dicom):
    # 查找图画的最大值和最小值
    vmin, vmax = tensor_dicom.min(), tensor_dicom.max()
    # 归一化
    tensor_dicom= (tensor_dicom- vmax ) / (max_val - vmin)
    
    return tensor_dicom

完成根据method句柄挑选预处理办法：

if method == "norm":
    # 归一化
    tensor_dcm = norm_tensor(tensor_dcm)
elif method == "uint8":
    # 归一化
    tensor_dcm = norm_tensor(tensor_dcm)
    # 量化
    tensor_dcm = (tensor_dcm * 255).clamp(0, 255).to(torch.uint8)

绘图

由于dicom图画为三维数据，可视化时咱们一般将其在z轴上分为多个切片顺次可视化，本文选用的办法是，选用5*5网格可视化至多25个切片。

def show_dciom(tensor_dicom):
    # 查找图画的最大最小值
    vmin, vmax = tensor_dicom.min(), tensor_dicom.max()
    
    # 创立一个图形窗口
    fig, axes = plt.subplots(5, 5, figsize=(15, 15))  # 5x5 网格布局

    count = 0
    length = tensor_dicom.size()[0]
    for i in range(25):
        if count < length:
            count += 1
        else:
            return
        # 获取当时图画的坐标
        ax = axes[i // 5, i % 5]
        # 显现图片
        ax.imshow(tensor_dicom[i], cmap='gray') # , vmin=vmin, vmax=vmax
        ax.axis('off')  # 封闭坐标轴
    
    plt.tight_layout() # 防止堆叠
    plt.title(f"Layer {i}")
    plt.show()

这里有一点需求比较留意，在ax.imshow()函数中，咱们指定了vmin和vmax参数；这是由于当该参数未被指守时，imshow函数将会主动调整点的亮度，使值最大的点对应255亮度，值最小的点对应0亮度。鉴于相邻切片最大、最小像素值或许存在较大差异，这将使得相邻切片的图画亮度较反常，如下图：

这两张图的左上角区域实际上亮度邻近，但从可视化图画来看，存在较大差异，这将对调查带来误解。

可视化频率散布直方图

可视化MRI图画的频率散布直方图在医学影像处理中有重要意义，首要包含以下几个方面：

• 图画对比度剖析：频率散布直方图能够显现MRI图画中不同灰度等级（或像素强度）的散布状况。经过剖析直方图的形状和规模，能够了解图画的对比度。例如，直方图的散布规模较广表明图画对比度较高，能够更好地区别不同安排或结构。

• 图画均衡化：经过直方图均衡化，能够改进图画的对比度，使得低对比度的区域愈加明晰。均衡化进程经过重新分配图画中的像素值，使得直方图的散布愈加均匀，然后增强图画的视觉作用。

• 安排切割：频率散布直方图能够协助确认恰当的阈值，以进行图画切割。经过剖析直方图，能够挑选适宜的阈值将不同安排或病变从布景中分离出来。

• 图画质量评价：直方图剖析能够提醒图画的质量问题，例如过暗或过亮的图画，或许图画噪声的影响。经过直方图的形状，能够评价图画是否需求进一步的处理或优化。

在制造频率散布直方图前，需求先将三维向量展平，本文选用plt.hist函数制造

def show_hist(tensor_dicom):
    # 将一切图片的像素值展平为一个一维数组
    pixel_values = tensor_dicom.numpy().flatten()

    # 制造直方图
    plt.figure(figsize=(10, 6))
    plt.hist(pixel_values, bins=50, color='gray', edgecolor='black')
    plt.title('Histogram of All Pixel Values')
    plt.xlabel('Pixel Value')
    plt.ylabel('Frequency')
    plt.grid(True)
    plt.show()

直方图呈现如下分步，在val=0邻近有一顶峰，这是由于MRI图画中大部分区域并不存在人体安排，为空值0。
假使除零以外的点过分会集在较小值（<100），那么很或许是由于MRI图画中呈现了一个亮度极大的噪点，使得以该噪点亮度为最值归一化质量较差，关于这种景象，能够用99%分位数替代最大值，并将99%分位数归一化至亮度为200. (比起归一化至255，这将答应亮度最大1%的像素点亮度值有区别）。
本例中图画质量均较高，故不需求做特别处理。

代码汇总

代码架构

主函数

# main.py
# Import custom utility functions
from utils import read_one_dicom, show_dciom, show_hist

# Define the directory containing the DICOM images
dicom_dir = "./train_images/4003253/1054713880"

# Read the DICOM image into a tensor with uint8 data type
tensor_dicom = read_one_dicom(dicom_dir, method="uint8")

# Display the DICOM image slices in a 5x5 grid layout
show_dciom(tensor_dicom)

# Plot the histogram of pixel values from the DICOM image slices
show_hist(tensor_dicom)

# Convert the tensor to a NumPy array for further processing or inspection
np_img = tensor_dicom.numpy()

包文件

from .preprocess import read_one_dicom

from .show import show_dciom
from .show import show_hist

读取&预处理

# preprocess.py
import numpy as np
import torch
import os
import re
import pydicom
from tqdm import tqdm

def norm_tensor(tensor_dicom):
    """
    Normalize the image tensor to the range [0, 1].

    Args:
        tensor_dicom (torch.Tensor): Tensor containing image data.

    Returns:
        torch.Tensor: Normalized image tensor.
    """
    # Calculate the maximum and minimum values of the image tensor
    vmin, vmax = tensor_dicom.min(), tensor_dicom.max()

    # Normalize the image tensor to the range [0, 1]
    tensor_dicom = (tensor_dicom - vmin) / (vmax - vmin)
    
    return tensor_dicom

def extract_number(filepath):
    """
    Extract the numeric part from the DICOM filename.

    Args:
        filepath (str): Path to the DICOM file.

    Returns:
        int: Extracted number from the filename. Returns float('inf') if not found.
    """
    # Get the filename (including extension)
    filename = os.path.basename(filepath)
    # Extract numeric part from filename, assuming filenames end with digits, e.g., '1.dcm'
    match = re.search(r'(\d+)\.dcm$', filename)
    return int(match.group(1)) if match else float('inf')

def read_one_dicom(dicom_dir, method = "", bar_title = ""):
    """
    Reads DICOM files from a directory and converts them into a PyTorch tensor.

    Args:
        dicom_dir (str): Directory containing DICOM files.
        method (str): Optional method to process the tensor ('norm' for normalization, 'uint8' for normalization and conversion to uint8).
        bar_title (str): Optional title for the progress bar.

    Returns:
        torch.Tensor: PyTorch tensor containing image data from DICOM files.
    """
    # Get all DICOM files and sort them based on numeric part of the filename
    dicom_files = [os.path.join(dicom_dir, f) for f in os.listdir(dicom_dir) if f.endswith('.dcm')]    
    dicom_files.sort(key=extract_number)

    # Create an empty list to store image data
    dcm_list = []

    # Initialize tqdm progress bar
    with tqdm(total=len(dicom_files), desc='Processing DICOM files', unit='dcm', unit_scale=True, unit_divisor=1000000) as pbar:
        # Iterate over each DICOM file and read image data
        for count, dcm_file in enumerate(dicom_files, start=1):
            # Read the DICOM file
            dcm = pydicom.dcmread(dcm_file)

            # Extract and convert image data to a NumPy array
            image_data = dcm.pixel_array.astype(np.float32)

            # Add the image data to the list
            dcm_list.append(image_data)

            # Update progress bar description
            pbar.set_description(bar_title + 'Reading')

            # Update progress bar
            pbar.update(1)

    # Convert the list of image data to a PyTorch tensor and stack into a 3D tensor
    tensor_dicom = torch.stack([torch.tensor(image_data) for image_data in dcm_list])

    if method == "norm":
        # Normalize the image tensor
        tensor_dicom = norm_tensor(tensor_dicom)
    elif method == "uint8":
        # Normalize the image tensor
        tensor_dicom = norm_tensor(tensor_dicom)
        # Scale the tensor values to the range [0, 255] and convert to uint8 type
        tensor_dicom = (tensor_dicom * 255).clamp(0, 255).to(torch.uint8)

    return tensor_dicom

可视化、制造直方图

# show.py
import numpy as np
import torch
import matplotlib.pyplot as plt

def show_dciom(tensor_dicom):
    """
    Display MRI image slices in a 5x5 grid layout.

    Parameters:
    tensor_dicom (torch.Tensor): Tensor containing MRI image slices, expected shape is (N, H, W),
                                 where N is the number of slices, and H and W are the height and width of the images.
    """
    # Calculate the minimum and maximum pixel values in the tensor
    vmin, vmax = tensor_dicom.min(), tensor_dicom.max()
    
    # Create a figure with a 5x5 grid layout
    fig, axes = plt.subplots(5, 5, figsize=(15, 15))  # 5x5 grid layout

    count = 0
    length = tensor_dicom.size(0)
    for i in range(25):
        if count < length:
            count += 1
        else:
            return
        # Get the current subplot's axis
        ax = axes[i // 5, i % 5]
        # Display the image
        ax.imshow(tensor_dicom[count - 1], cmap='gray', vmin=vmin, vmax=vmax)
        ax.axis('off')  # Hide the axis
    
    plt.tight_layout()  # Adjust layout to prevent overlap
    plt.title(f"Layer {i + 1}")  # Title indicating the last displayed slice
    plt.show()

def show_hist(tensor_dicom):
    """
    Plot the histogram of pixel values for all MRI image slices.

    Parameters:
    tensor_dicom (torch.Tensor): Tensor containing MRI image slices, expected shape is (N, H, W).
    """
    # Flatten all image pixel values into a single 1D array
    pixel_values = tensor_dicom.numpy().flatten()

    # Plot the histogram
    plt.figure(figsize=(10, 6))
    plt.hist(pixel_values, bins=50, color='gray', edgecolor='black')
    plt.title('Histogram of All Pixel Values')
    plt.xlabel('Pixel Value')
    plt.ylabel('Frequency')
    plt.grid(True)
    plt.show()

下篇预告

评论本题的解题办法

制造不易，请帮我点一个免费的赞，谢谢！