【图像去噪】基于傅里叶的稳态去MRI共振伪影matlab1仿真

✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。

🍎 往期回顾关注个人主页：Matlab科研工作室

🍊个人信条：格物致知,完整Matlab代码及仿真咨询内容私信。

🔥 内容介绍

磁共振成像（MRI）凭借无辐射、软组织分辨率高的优势，已成为神经系统、肌肉骨骼系统等疾病诊断的核心影像技术。然而，MRI 成像过程中易受多种因素干扰产生伪影，其中共振伪影（又称周期性伪影）是最常见且影响诊断的类型之一。其本质是 MRI 信号采集时，射频脉冲激励、梯度场切换或患者生理运动（如呼吸、心跳）导致的信号周期性波动，在图像上表现为沿特定方向（多为频率编码方向）分布的明暗交替条纹，严重时会掩盖病灶细节 —— 例如，脑部 MRI 中，共振伪影可能与脑白质病变信号混淆，导致误诊率升高 15%-20%；在膝关节 MRI 中，伪影可能掩盖半月板撕裂的细微结构，影响手术方案制定。

从成因来看，MRI 共振伪影主要分为两类：

1. 设备源性共振伪影：射频线圈阻抗不匹配、梯度场非线性畸变，导致信号在频率域出现周期性噪声峰值，映射到图像域形成条纹；

1. 生理源性共振伪影：患者呼吸、心跳等周期性运动，使组织位置随时间周期性变化，导致采集的 K 空间数据出现周期性偏差，最终在图像上呈现与运动频率一致的伪影。

传统 MRI 去噪方法（如高斯滤波、中值滤波）在处理共振伪影时存在明显局限：

• 空间域滤波的盲目性：高斯滤波虽能平滑噪声，但会模糊病灶边缘（如肿瘤边界），中值滤波对周期性条纹伪影抑制效果有限，且易破坏图像纹理细节；

• 缺乏频率针对性：共振伪影在频率域具有明确的周期性特征（表现为离散的噪声峰值），而传统方法未利用这一特性，去噪效率低且易引入新的图像失真；

• 动态适应性差：不同部位、不同扫描参数下的 MRI 共振伪影频率、强度差异大，传统固定参数的滤波方法难以实现稳态去噪，需反复调整参数，适用性受限。

傅里叶变换（Fourier Transform, FT）作为连接 “空间域” 与 “频率域” 的核心工具，为 MRI 共振伪影的稳态去除提供了理想解决方案。其核心优势在于：能将空间域的周期性伪影转化为频率域的离散噪声峰值，通过针对性的频率域滤波实现伪影精准去除，同时最大程度保留有效图像信号，实现 “精准去噪 - 细节保留” 的稳态平衡。

⛳️ 运行结果

📣 部分代码

(model.protonDensity))

    model.size([2 1 3]) = size3(model.protonDensity); % XYZ

elseif (~isscalar(model.deltaB0))

    model.size([2 1 3]) = size3(model.deltaB0); % XYZ

elseif (~isscalar(model.T2))

    model.size([2 1 3]) = size3(model.T2); % XYZ

else

    error('forecast: Model is all scalars, at least one of the fields protonDensity, deltaB0, or T2 must be a matrix');

end

% Check that if protonDensity/deltaB0/T2 are given as matrices that the sizes match

if (~isscalar(model.deltaB0) && any(size3(model.deltaB0) ~= model.size([2 1 3])))

    error('forecast: protonDensity and deltaB0 are given as matrix, but their sizes do not match');

end

if (~isscalar(model.T2) && any(size3(model.T2) ~= model.size([2 1 3])))

    error('forecast: protonDensity and/or deltaB0 and T2 are given as matrix, but their sizes do not match');

end

if (~isfield(model, 'resolution'))

    error('forecast: Model structure is missing the resolution field');

end

model.fieldOfView = model.size .* model.resolution; % XYZ

%% Acquisition parameters

if (~isstruct(acquisition))

    error('Acquisition parameter must be a structure')

end

if (~isfield(acquisition, 'kspaceSamplingTimes'))

    error('forecast: Acquisition structure is missing the kspaceSamplingTimes field');

end

if (~isfield(acquisition, 'kspaceSamplingTimesRefocused'))

    acquisition.kspaceSamplingTimesRefocused = acquisition.kspaceSamplingTimes;

end

if (~isequal(size(acquisition.kspaceSamplingTimes), size(acquisition.kspaceSamplingTimesRefocused)))

    error('forecast: The size of kspaceSamplingTimes must be equal to the size kspaceSamplingTimesRefocused');

end

if (~isfield(acquisition, 'resolution'))

    error('forecast: Acquisition structure is missing the resolution field');

end

if (any(model.resolution > acquisition.resolution))

    error('forecast: Model resolution must be higher than acquisition resolution')

end

acquisition.size([2 1 3]) = size3(acquisition.kspaceSamplingTimes); % XYZ

acquisition.fieldOfView = acquisition.size .* acquisition.resolution; % XYZ

%% Options parameters

if (~isstruct(options))

    error('Options parameter must be a structure');

end

if (~isfield(options, 'verbose'))

    options.verbose = false;

end

if (~isfield(options, 'forceGeneral'))

    options.forceGeneral = false;

end

% Set gyromagnetic ratio constant

gm = 42.57747892e6; % Hz / T

kspaceData = zeros(acquisition.size([2 1 3]));

for I=1:size(Ts,1)

    t_B0 = Ts(I,1);

    t_T2 = Ts(I,2);

    if (options.verbose)

        fprintf('Simulating time %d / %d: %f ms / %f ms\n', I, size(Ts,1), t_B0, t_T2);

    end

    % Analytical description of transverse magnetization without encoding at time t

    Mtrans = model.protonDensity .* exp(-2i * pi * gm * model.deltaB0 * (t_B0/1000) - t_T2 ./ model.T2);

    % Spatially encode in all dimensions

    kSimulation = fftc(Mtrans);

    % Crop k-space to match scan resolution

    k = crop(kSimulation, round(model.size([2 1 3]) .* scale([2 1 3])));

    % Correct intensity

    k = k / prod(sqrt(size3(kSimulation) ./ round(model.size([2 1 3]) .* scale([2 1 3]))));

    % Correct field of view by either cropping or zero-padding in image-space

    k = fftc(zeroPadOrCrop(ifftc(k), acquisition.size([2 1 3])));

    % Select elements from k-space that were encoded at time t and save

    % them in the final simulated k-space

    kspaceData(Tsind==I) = k(Tsind==I);

end

end

% Size function that ensures the returned number of dimensions is at least 3

function D = size3 (M)

D = size(M);

if (length(D) < 3)

    D(3) = 1;

end

end

🔗 参考文献

[1]赵敬聪.屈光参差性弱视儿童视觉运动功能损害的功能磁共振成像研究[D].天津医科大学,2011.DOI:10.7666/d.y2001035.

🎈 部分理论引用网络文献，若有侵权联系博主删除

 👇 关注我领取海量matlab电子书和数学建模资料 

🏆团队擅长辅导定制多种科研领域MATLAB仿真，助力科研梦：

🌟 各类智能优化算法改进及应用

生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化、背包问题、 风电场布局、时隙分配优化、 最佳分布式发电单元分配、多阶段管道维修、 工厂-中心-需求点三级选址问题、 应急生活物质配送中心选址、 基站选址、 道路灯柱布置、 枢纽节点部署、 输电线路台风监测装置、 集装箱调度、 机组优化、 投资优化组合、云服务器组合优化、 天线线性阵列分布优化、CVRP问题、VRPPD问题、多中心VRP问题、多层网络的VRP问题、多中心多车型的VRP问题、 动态VRP问题、双层车辆路径规划（2E-VRP）、充电车辆路径规划（EVRP）、油电混合车辆路径规划、混合流水车间问题、 订单拆分调度问题、 公交车的调度排班优化问题、航班摆渡车辆调度问题、选址路径规划问题、港口调度、港口岸桥调度、停机位分配、机场航班调度、泄漏源定位、冷链、时间窗、多车场等、选址优化、港口岸桥调度优化、交通阻抗、重分配、停机位分配、机场航班调度、通信上传下载分配优化

🌟 机器学习和深度学习时序、回归、分类、聚类和降维

2.1 bp时序、回归预测和分类

2.2 ENS声神经网络时序、回归预测和分类

2.3 SVM/CNN-SVM/LSSVM/RVM支持向量机系列时序、回归预测和分类

2.4 CNN|TCN|GCN卷积神经网络系列时序、回归预测和分类

2.5 ELM/KELM/RELM/DELM极限学习机系列时序、回归预测和分类

2.6 GRU/Bi-GRU/CNN-GRU/CNN-BiGRU门控神经网络时序、回归预测和分类

2.7 ELMAN递归神经网络时序、回归\预测和分类

2.8 LSTM/BiLSTM/CNN-LSTM/CNN-BiLSTM/长短记忆神经网络系列时序、回归预测和分类

2.9 RBF径向基神经网络时序、回归预测和分类

2.10 DBN深度置信网络时序、回归预测和分类

2.11 FNN模糊神经网络时序、回归预测

2.12 RF随机森林时序、回归预测和分类

2.13 BLS宽度学习时序、回归预测和分类

2.14 PNN脉冲神经网络分类

2.15 模糊小波神经网络预测和分类

2.16 时序、回归预测和分类

2.17 时序、回归预测预测和分类

2.18 XGBOOST集成学习时序、回归预测预测和分类

2.19 Transform各类组合时序、回归预测预测和分类

方向涵盖风电预测、光伏预测、电池寿命预测、辐射源识别、交通流预测、负荷预测、股价预测、PM2.5浓度预测、电池健康状态预测、用电量预测、水体光学参数反演、NLOS信号识别、地铁停车精准预测、变压器故障诊断

🌟图像处理方面

图像识别、图像分割、图像检测、图像隐藏、图像配准、图像拼接、图像融合、图像增强、图像压缩感知

🌟 路径规划方面

旅行商问题（TSP）、车辆路径问题（VRP、MVRP、CVRP、VRPTW等）、无人机三维路径规划、无人机协同、无人机编队、机器人路径规划、栅格地图路径规划、多式联运运输问题、 充电车辆路径规划（EVRP）、 双层车辆路径规划（2E-VRP）、 油电混合车辆路径规划、 船舶航迹规划、 全路径规划规划、 仓储巡逻、公交车时间调度、水库调度优化、多式联运优化

🌟 无人机应用方面

无人机路径规划、无人机控制、无人机编队、无人机协同、无人机任务分配、无人机安全通信轨迹在线优化、车辆协同无人机路径规划、

🌟 通信方面

传感器部署优化、通信协议优化、路由优化、目标定位优化、Dv-Hop定位优化、Leach协议优化、WSN覆盖优化、组播优化、RSSI定位优化、水声通信、通信上传下载分配

🌟 信号处理方面

信号识别、信号加密、信号去噪、信号增强、雷达信号处理、信号水印嵌入提取、肌电信号、脑电信号、信号配时优化、心电信号、DOA估计、编码译码、变分模态分解、管道泄漏、滤波器、数字信号处理+传输+分析+去噪、数字信号调制、误码率、信号估计、DTMF、信号检测

🌟电力系统方面

微电网优化、无功优化、配电网重构、储能配置、有序充电、MPPT优化、家庭用电、电/冷/热负荷预测、电力设备故障诊断、电池管理系统（BMS）SOC/SOH估算（粒子滤波/卡尔曼滤波）、 多目标优化在电力系统调度中的应用、光伏MPPT控制算法改进（扰动观察法/电导增量法）、电动汽车充放电优化、微电网日前日内优化、储能优化、家庭用电优化、供应链优化

🌟 元胞自动机方面

交通流 人群疏散 病毒扩散 晶体生长 金属腐蚀

🌟 雷达方面

卡尔曼滤波跟踪、航迹关联、航迹融合、SOC估计、阵列优化、NLOS识别

🌟 车间调度

零等待流水车间调度问题NWFSP 、 置换流水车间调度问题PFSP、 混合流水车间调度问题HFSP 、零空闲流水车间调度问题NIFSP、分布式置换流水车间调度问题 DPFSP、阻塞流水车间调度问题BFSP

👇

5 往期回顾扫扫下方二维码