3D MRI运动伪影模拟

3D MRI运动伪影模拟

背景与动机

在医学影像中，运动伪影是一种常见的问题，特别是在MRI图像中。这种伪影会影响医生的诊断准确性。为了解决这一挑战，我们开发了一个伪影模拟工具，帮助研究人员分析和减轻伪影的影响。

核心功能

该项目的核心是通过随机模拟运动轨迹，结合FFT变换，在频域中生成运动伪影。

代码实现了3D MRI图像伪影的模拟，涉及图像频域变换、随机运动轨迹的模拟以及伪影生成的逻辑。是一个医学成像中的运动伪影生成项目。核心功能集中在利用FFT、随机运动模拟和空间变换算法实现图像的随机扰动。

关键因素：

1. 代码主要用于MRI图像的运动伪影模拟。
2. 涉及到随机生成的平移和旋转运动参数，并通过FFT变换在频域中引入伪影。
3. 数学原理以及如何从频域变换到时域。

实现算法

1. 图像处理基础：FFT和iFFT

通过FFT（快速傅里叶变换）将空间域图像转换到频域。

• 数学表达式：
  $$F(u,v,w)=\sum _{x=0}^{N-1}\sum _{y=0}^{N-1}\sum _{z=0}^{N-1}f(x,y,z)e^{-j2\pi (ux/N+vy/N+wz/N)}$$

  • f(x, y, z) 是输入图像
  • F(u, v, w) 是频域表示

2. 运动伪影的随机轨迹模拟

随机生成平移和旋转参数：

• 平移模拟：利用正态分布生成随机平移值，并通过相位变化引入伪影。
• 旋转模拟：基于3D旋转矩阵构造旋转效果。

3. 中心k-space保护机制

在医学成像中，k-space的中心包含重要信息。因此，代码通过get_center_rect函数保护中心区域，防止过度伪影。

4. 频域中的随机扰动

通过_translate_freq_domain函数，将平移转换为频域中的相位变化：

• 数学推导：
  $$F^{\prime }(u,v,w)=F(u,v,w)\cdot e^{-j2\pi ({\Delta }_{x}u+{\Delta }_{y}v+{\Delta }_{z}w)}$$

5. 伪影合成

结合FFT与逆FFT，将频域中的伪影引入时域图像。

代码结构解析

核心模块

1. MotionSimLayer类
   • 参数说明：
     • std_rotation_angle：旋转角度的标准差
     • std_translation：平移的标准差
     • corrupt_pct_range：伪影覆盖百分比范围
   • 功能：
     模拟随机运动轨迹，并将伪影引入图像。
2. 辅助函数
   • create_rotation_matrix_3d：生成3D旋转矩阵
   • segment_array_by_locs：生成随机分段
   • _fft_im与_ifft_im：FFT与iFFT转换

三维伪影模拟的创新点及突破

1. 三维与二维方法的对比

二维伪影模拟通常基于平面操作，考虑水平和垂直方向的运动对图像的影响。然而，MRI成像是三维的，伪影不仅存在于单一切片，还会在深度方向上影响图像的整体质量。因此，二维方法在以下方面存在局限性：

• 运动复杂性：无法模拟深度方向（z轴）的运动伪影。
• 真实感不足：二维方法只反映平面伪影，无法重现MRI中的全局伪影特性。
• 局部特性偏重：无法模拟跨切片的连续性伪影。

2. 三维方法的创新与突破

本算法在以下几个方面突破了二维伪影的局限性：

（1）三维旋转与平移模拟

核心创新：引入了三维旋转矩阵与平移矢量，全面模拟三维空间中的运动。

• 三维旋转矩阵生成：
  通过create_rotation_matrix_3d函数，根据随机角度生成三个方向上的旋转矩阵，并组合为最终的三维旋转矩阵：
  $$R=R_x\cdot R_y\cdot R_z$$
  其中：
  $$R_x=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\theta }_x& \sin {\theta }_x\\ 0& -\sin {\theta }_x& \cos {\theta }_x\end{array}\right],R_y=\left[\begin{array}{ccc}\cos {\theta }_y& 0& -\sin {\theta }_y\\ 0& 1& 0\\ \sin {\theta }_y& 0& \cos {\theta }_y\end{array}\right],R_z=\left[\begin{array}{ccc}\cos {\theta }_z& \sin {\theta }_z& 0\\ -\sin {\theta }_z& \cos {\theta }_z& 0\\ 0& 0& 1\end{array}\right]$$
  这种方法显著提升了对三维运动伪影的模拟能力。

• 三维平移模拟：在频域中引入与三维平移对应的相位偏移：
  $$F^{\prime }(u,v,w)=F(u,v,w)\cdot e^{-j2\pi ({\Delta }_{x}u+{\Delta }_{y}v+{\Delta }_{z}w)}$$
  通过独立控制x、y、z三个方向的平移幅度，能够生成更加复杂的伪影效果。

（2）三维频域变换与伪影生成

通过3D FFT与iFFT，支持从三维频域到空间域的转换：

• 三维FFT公式：
  $$F(u,v,w)=\sum _{x=0}^{N-1}\sum _{y=0}^{M-1}\sum _{z=0}^{L-1}f(x,y,z)e^{-j2\pi (ux/N+vy/M+wz/L)}$$

  • FFT用于将三维图像变换到频域，方便引入运动扰动。
  • iFFT将频域图像转换回空间域，重建伪影图像。

• 保护中心k-space：为避免过度伪影对图像信息的完全破坏，设计了三维中心保护机制（get_center_rect函数扩展到三维）。

（3）引入非均匀FFT（NUFFT）

在复杂旋转与空间变换中，传统FFT可能无法高效处理非均匀采样点。本算法支持NUFFT，能够处理不规则运动轨迹，并提升伪影模拟的灵活性。但调试失败，没有实现

（4）随机轨迹与伪影分布控制

• 分段伪影模拟：通过segment_array_by_locs和assign_segments_to_random_indices，实现对伪影的随机化分布与动态模拟。
• 伪影强度可控：通过corrupt_pct_range参数，灵活调整伪影的覆盖范围，从局部影响到全局影响均可模拟。

项目应用与实现

应用场景

本项目的主要应用领域是医学影像模拟，特别是MRI图像中的运动伪影研究。它适用于：

1. 算法验证：开发和测试运动伪影校正算法。
2. 图像分析：研究运动伪影对图像质量和诊断准确性的影响。
3. 教育与培训：为医学影像学的教育与研究提供真实的运动伪影数据。

项目实现流程

该项目采用以下步骤完成影像伪影的生成和处理：

1. 加载DICOM影像数据

通过load_dicom_as_3d函数，将DICOM格式的切片加载并堆叠为三维数组。此步骤对切片进行排序并提取像素间距信息，确保影像的空间一致性。

2. 三维伪影生成

核心算法由MotionSimLayer类实现：

• 在加载的三维图像上模拟随机的运动伪影。
• 通过FFT和随机运动轨迹生成伪影，模拟出更接近真实场景的三维MRI影像。

3. 保存处理结果

将处理后的三维影像保存为新的DICOM切片，支持与其他影像处理系统兼容的后续使用。

项目依赖

本项目基于Ben Duffy的motion-correction项目，进行了以下扩展和优化：

1. 增强了伪影生成逻辑，支持更复杂的运动模式和三维伪影特性。
2. 提供了一个完整的工作流，从数据加载到保存，涵盖伪影模拟的全流程。
3. 引入了灵活的配置参数（如std_rotation_angle和corrupt_pct_range），适应多种研究需求。

import os
import numpy as np
import pydicom
from clear import MotionSimLayer

# === Step 1: 加载 DICOM 文件并堆叠为 3D ===
def load_dicom_as_3d(directory):
    """从目录加载 DICOM 切片并堆叠为 3D 数组"""
    dicom_files = [os.path.join(directory, f) for f in os.listdir(directory) if f.endswith('.dcm')]
    slices = [pydicom.dcmread(f) for f in dicom_files]
    slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))  # 按切片位置排序

    pixel_arrays = [s.pixel_array for s in slices]
    voxel_spacing = (
        float(slices[0].PixelSpacing[0]),
        float(slices[0].PixelSpacing[1]),
        float(abs(slices[1].ImagePositionPatient[2] - slices[0].ImagePositionPatient[2]))
    )

    return np.stack(pixel_arrays, axis=-1), slices, voxel_spacing

# === Step 2: 使用 MotionSimLayer 处理 3D 数据 ===
def process_3d_with_motion_sim(input_3d, motion_sim_layer):
    """对输入的 3D 数组应用 MotionSimLayer"""
    processed_3d = motion_sim_layer(input_3d)
    return processed_3d

# === Step 3: 将处理后的数据保存为 DICOM 和 3D 格式 ===
def save_processed_as_dicom(processed_3d, original_slices, output_dir):
    """将处理后的 3D 数组保存为 DICOM 切片"""
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    for i in range(processed_3d.shape[2]):
        slice_2d = processed_3d[:, :, i].astype(np.int16)
        dicom_slice = original_slices[i]
        dicom_slice.PixelData = slice_2d.tobytes()
        dicom_slice.save_as(os.path.join(output_dir, f"slice_{i:04d}.dcm"))




# === Main 执行部分 ===
if __name__ == "__main__":
    # 输入参数
    dicom_input_dir = "./dicom_input"  # DICOM 输入目录
    dicom_output_dir = "./dicom_output"  # DICOM 输出目录


    # 加载 DICOM 数据
    print("Loading DICOM files...")
    input_3d, original_slices, voxel_spacing = load_dicom_as_3d(dicom_input_dir)
    print(f"Loaded DICOM data with shape: {input_3d.shape}")

    # 初始化 MotionSimLayer 类
    motion_sim_layer = MotionSimLayer(
        std_rotation_angle=5,  # 旋转标准差，控制随机旋转的幅度
        std_translation=10,  # 平移标准差，控制随机平移的幅度
        corrupt_pct_range=(10, 20),  # 伪影影响范围的百分比（范围内随机选择）
        freq_encoding_dim=(0, 1, 2),  # 随机选择频率编码的维度
        preserve_center_pct=0.07,  # k空间中心保护区域的百分比
        apply_mask=True,  # 是否保护中心频率
        nufft=False,  # FFT（NUFFT），调试失败
        corruption_scheme='piecewise_constant',  # 伪影模式，可选值有 'piecewise_constant', 'piecewise_transient', 'gaussian'
        n_seg=8,  # 最大分段数量（用于分段伪影模式）
        fixed_n_seg=False,  # 是否使用固定分段数量

    )


    print("Processing 3D data with MotionSimLayer...")
    processed_3d = process_3d_with_motion_sim(input_3d, motion_sim_layer)
    print("Processing complete.")

    # 保存为 DICOM 切片
    print(f"Saving processed data as DICOM to {dicom_output_dir}...")
    save_processed_as_dicom(processed_3d, original_slices, dicom_output_dir)