【MRI】SENSE (Sensitivity Encoding) 算法仿真实验与原理剖析 (Matlab 实现)

何处闻韶

已于 2022-05-11 16:52:39 修改

阅读量5k

点赞数 24

分类专栏：【智慧医疗】文章标签： matlab 算法机器学习

于 2020-11-04 09:58:40 首次发布

本文链接：https://blog.csdn.net/qq_39478403/article/details/109278019

版权

【智慧医疗】专栏收录该内容

2 篇文章 10 订阅

订阅专栏

1. 加载全采样 MR 及其敏感度图并显示

1.1 预扫描数据显示

1.2 根据预扫描数据初次生成敏感度图并显示

2. 获取部分参数

3. 将全采样 MR 转换为全采样 k 空间

4. 显示各线圈及 RSOS 组合重建的总体全采样 k 空间频谱

4.1 正常步骤

4.2 根据预扫描数据的 IFFT 再次生成敏感度图并显示

5. 设置加速因子

6. 构造欠采样等距掩模 (equispaced mask) 并显示

7. 模拟/生成欠采样 k 空间数据并显示

8. 欠采样 k 空间数据 IFFT 转换到图像域并显示

功能函数 rsos.m (root-sum-of-square)

function image = rsos(data)
% RSOS Root Sum of Squares Function.
% 
% The root Sum of Squares function necessary for converting 3D multichannel 
% complex data into 2D real valued data.
assert(length(size(data)) == 2 || length(size(data)) == 3, 'Data must be either 2D or 3D.');
image = abs(data) .^ 2;
image = sum(image, 3);
image = sqrt(image);
assert(length(size(image)) == 2);
end

%% Assignment 2-1: SENSE Reconstruction
%% SENSE
% (SENSitivity Encoding) parallel imaging reconstruction method.
% 
% See <http://mriquestions.com/senseasset.html here> for a basic introduction. 
% See <https://www.researchgate.net/publication/8178233_SMASH_SENSE_PILS_GRAPPA_how_to_choose_the_optimal_method 
% this review paper> (p.226) for a more detailed explanation.
% 
% The notation used in the code below is based on the review paper.

约定：相位编码方向 k_y 自底向上 (↑)、读出方向/频率编码方向 k_x 从左往右 (→)。

1. 加载全采样 MR 及其敏感度图并显示

1.1 预扫描数据显示

brainCoils 是 5 通道的大脑 MR 图像，分别来自 Nc = 5 个线圈各自的全采样结果，每张图像尺寸为 120 ×128。

coilMap 是 5 通道的线圈敏感度图，与 brainCoils 的 shape 或 size 相同，元素值一一对应。该敏感度图来自 预扫描 过程，(本实现暂未涉及)，默认正常工作时每次扫描环境下敏感度图相同 (本实现仅用于展示说明，降低严谨性是可容忍的)。

%% Loading data
brainCoilsData = load('brain_coil.mat');
brainCoils = brainCoilsData.brain_coil_tmp;

coilMapData = load('coil_sensitivity_map.mat');
coilMap = coilMapData.coil_map;  % 120×128×5 double

接着，分别显示各线圈的全采样 MR 图 brainCoils(:, :, Nc)：

%% Examine full sampled MR images of each coil
figure;
colormap('gray');

subplot(2, 3, 1);
imagesc(brainCoils(:,:,1));
title('coil-1 full sampled MR image');
colorbar();

subplot(2, 3, 2);
imagesc(brainCoils(:,:,2));
title('coil-2 full sampled MR image');
colorbar();

subplot(2, 3, 3);
imagesc(brainCoils(:,:,3));
title('coil-3 full sampled MR image');
colorbar();

subplot(2, 3, 4);
imagesc(brainCoils(:,:,4));
title('coil-4 full sampled MR image');
colorbar();

subplot(2, 3, 5);
imagesc(brainCoils(:,:,5));
title('coil-5 full sampled MR image');
colorbar();

然后，分别显示各线圈的敏感度图 coilMap(:, :, Nc)：

%% Examine sensitivity map of each coil
figure;
colormap('gray');

subplot(2, 3, 1);
imagesc(coilMap(:,:,1));
title('coil-1 sensitivity map');
colorbar();

subplot(2, 3, 2);
imagesc(coilMap(:,:,2));
title('coil-2 sensitivity map');
colorbar();

subplot(2, 3, 3);
imagesc(coilMap(:,:,3));
title('coil-3 sensitivity map');
colorbar();

subplot(2, 3, 4);
imagesc(coilMap(:,:,4));
title('coil-4 sensitivity map');
colorbar();

subplot(2, 3, 5);
imagesc(coilMap(:,:,5));
title('coil-5 sensitivity map');
colorbar();

1.2 根据预扫描数据初次生成敏感度图并显示

%% test about generating sensitivity map manually
BodybrainCoils = rsos(brainCoils);  % RSOS reconstruction for body MR image

figure;
colormap('gray');

subplot(2, 3, 1);
sensitivity_map_1 = brainCoils(:,:,1) ./ BodybrainCoils;
imagesc(sensitivity_map_1);
title('coil-1 sensitivity map (manual)');
colorbar();

subplot(2, 3, 2);
sensitivity_map_2 = brainCoils(:,:,2) ./ BodybrainCoils;
imagesc(sensitivity_map_2);
title('coil-2 sensitivity map (manual)');
colorbar();

subplot(2, 3, 3);
sensitivity_map_3 = brainCoils(:,:,3) ./ BodybrainCoils;
imagesc(sensitivity_map_3);
title('coil-3 sensitivity map (manual)');
colorbar();

subplot(2, 3, 4);
sensitivity_map_4 = brainCoils(:,:,4) ./ BodybrainCoils;
imagesc(sensitivity_map_4);
title('coil-4 sensitivity map (manual)');
colorbar();

subplot(2, 3, 5);
sensitivity_map_5 = brainCoils(:,:,5) ./ BodybrainCoils;
imagesc(sensitivity_map_5);
title('coil-5 sensitivity map (manual)');
colorbar();

subplot(2, 3, 6);
imagesc(BodybrainCoils);
title('RSOS reconstruction body MR image');
colorbar();

对比可见，根据公式手动生成的敏感度图不仅具有一些噪点瑕疵，而且对比度也存在差异。虽然 1.1 中预扫描给出的敏感度图可能是经过后处理的，但根本原因在于这种计算方式不符合实际，毕竟 MRI 系统实际先采集到的应该是 raw k-space 而非 MR image。为此，真正的计算方式详见 4.2 节。

敏感度图的计算

计算线圈敏感度是 SENSE 过程中初始而最重要的一步。低分辨率图像是在 全 FOV 下从各表面线圈分别采集的。这些表面线圈图像通过将其除以低分辨率总体线圈图像进行归一化。然后对数据进行 滤波、阈值化和点估计，以生成线圈敏感度图（如右图所示）。这些映射量化每个线圈接收区域内不同原点的信号的相对权重。

一旦线圈敏感度图被计算出来，正式工作的磁共振脉冲序列就开始了。

注意，SENSE 假定后续每次扫描时，敏感度图都是不变的。然而，事实上每次扫描的系统环境可能都有所不同，敏感度图会随扫描次数而发生 (细节性的) 变化。

2. 获取部分参数

phaseLength 表示 k 空间相位编码方向长度 ky_FOV = 120。(图像域 MR 图像/视野高度 FOV_y = 120)

freqLength 表示 k 空间频率编码方向长度 kx_FOV = 128。(图像域 MR 图像/视野宽度 FOV_x = 128)

coilNum 表示采用线圈数 Nc = 5。

%% size
[phaseLength, freqLength, coilNum] = size(brainCoils);  % ky_FOV=120, kx_FOV=128, Nc=5
disp(size(brainCoils));  % Phases, Frequencies, Coil Number

3. 将全采样 MR 转换为全采样 k 空间

fullKpace 为各线圈全采样 MR 图像对应的 k 空间域数据，size = 120×128×5，type = complex double

%% Convert to Fourier Domain with 2D Fast Fourier Transform
% Remember to use fftshift and ifftshift when using the Fourier Transform. 
% Otherwise, the k-space will not be centered properly.

fullKspace = ifftshift(fft2(fftshift(brainCoils)));  % 注意转换到 k 空间时, FFT 前后要分别 shift

fftshift 与 ifftshift

注意，当空域中的点 f(x,y) 移位时，在频域中只发生相移，并不影响其傅里叶变换的幅度。反之，当频域中的点 F(u,v) 移位时，相应的 f(x,y) 在空域中也只发生相移，不产生幅值变化。根据平移性质，为更清楚查看频率域的 k 空间数据的频谱 (复数值数据的模)，在把画面分成四分的基础上，可以进行换位/移位 (shift)，从而使频域原点 (直流成分) 平移到频谱中央。如下所示：

举个例子：

4. 显示各线圈及 RSOS 组合重建的总体全采样 k 空间频谱

4.1 正常步骤

fullKspaceImage 是通过 平方和根 (root-sum-of-square) 方法 组合所有线圈各自的全采样 k 空间 复数值数据 重建的总体全采样 k 空间 实数值频谱，size = 120×128，type = double

rsos(fullKspace(:, :, Nc) 则分别为由各线圈各自的全采样 k 空间 复数值数据 取模得到的 实数值频谱，size = 120×128，type = double

%% Examine raw k-space.
fullKspaceImage = rsos(fullKspace);  % RSOS 重建得到实值全采样 k 空间频谱图像

figure;
colormap('gray');

subplot(2, 3, 1);
imagesc(rsos(fullKspace(:,:,1)), [0, 10000]);
title('coil-1 k-space');
colorbar();

subplot(2, 3, 2);
imagesc(rsos(fullKspace(:,:,2)), [0, 10000]);
title('coil-2 k-space');
colorbar();

subplot(2, 3, 3);
imagesc(rsos(fullKspace(:,:,3)), [0, 10000]);
title('coil-3 k-space');
colorbar();

subplot(2, 3, 4);
imagesc(rsos(fullKspace(:,:,4)), [0, 10000]);
title('coil-4 k-space');
colorbar();

subplot(2, 3, 5);
imagesc(rsos(fullKspace(:,:,5)), [0, 10000]);
title('coil-5 k-space');
colorbar();

subplot(2, 3, 6);
imagesc(fullKspaceImage, [0, 10000]);
title('RSOS reconstruction k-space');
colorbar();

4.2 根据预扫描数据的 IFFT 再次生成敏感度图并显示

事实上，MRI 系统采集到的本应就是 raw k-space 数据，只不过代码给的测试文件是图像域 MR 图像而已。

%% =================================== 专门测试敏感度图 =====================================
rawfullKspace = fullKspace;
acsBrainCoils = ifftshift(ifft2(fftshift(rawfullKspace)));  % 128×120×5 complex double
acsImage = rsos(acsBrainCoils);  % RSOS 重建 128×120 double

figure;
colormap('gray');

subplot(2, 3, 1);
sens_map_1 = rsos(acsBrainCoils(:,:,1)) ./ acsImage;
imagesc(sens_map_1);
title('coil-1 sensitivity map (manual)');
colorbar();

subplot(2, 3, 2);
sens_map_2 = rsos(acsBrainCoils(:,:,2)) ./ acsImage;
imagesc(sens_map_2);
title('coil-2 sensitivity map (manual)');
colorbar();

subplot(2, 3, 3);
sens_map_3 = rsos(acsBrainCoils(:,:,3)) ./ acsImage;
imagesc(sens_map_3);
title('coil-3 sensitivity map (manual)');
colorbar();

subplot(2, 3, 4);
sens_map_4 = rsos(acsBrainCoils(:,:,4)) ./ acsImage;
imagesc(sens_map_4);
title('coil-4 sensitivity map (manual)');
colorbar();

subplot(2, 3, 5);
sens_map_5 = rsos(acsBrainCoils(:,:,5)) ./ acsImage;
imagesc(sens_map_5);
title('coil-5 sensitivity map (manual)');
colorbar();

subplot(2, 3, 6);
imagesc(acsImage);
title('RSOS reconstruction body MR image');
colorbar();

可以看到，本次生成的敏感度图与程序初始输入的敏感度图十分接近了！

5. 设置加速因子

downSamplingRate 为下/欠采样率，又称缩减因子 (reduction factor) / 加速因子 (accelation factor)，符号通常为 R。此处设置 R = 5 表示将进行 5 倍加速 MRI，这在笛卡尔 k 空间顺序采集中，意味着每 R = 5 条相位编码线的位置才真正地采集到 1 条 (acquired line)，其位置将由 1 值 mask 覆盖表示已采集；其余 4 条 (unacquired line) 的位置将由 0 值 mask 覆盖表示未采集。从而，实现 R = 5 倍欠采样 / 加速。

%% setting reduction/accelation factor
downSamplingRate = 5;  % 缩减/加速因子 R = 5 —— 5 倍欠采样/加速

6. 构造欠采样等距掩模 (equispaced mask) 并显示

fullKpace 为各线圈全采样 MR 图像对应的 k 空间域数据，size = 120×128×5，type = complex double

mask 为与 fullKpace 等 size 和 type 的 0/1 等距掩模 (不是随机掩模)，用于模拟/制作 R = 5 倍欠采样 k 空间数据。

%% Making the mask
% Creating a mask for the downsampling and the ACS lines.

mask = zeros(size(fullKspace));  % Making a mask for the brain coils.
for line = 1:phaseLength   % Goes along the phase encoding axis. (vertical direction ky - dim=0)
    if rem(line, downSamplingRate) == 0  % r = rem(a,b) returns the remainder after division of a by b, 
        mask(line, :, :) = 1;  % Broadcasting the value of 1 to mask(i, :, :).  % 每 R 行采集一条相位编码线, mask=1
    end    
end
assert(isequal(size(mask), size(brainCoils)), 'Mask size is incorrect.')  % mask 要和 full FOV MR image 等尺寸

然后，显示构造好的 R = 5 倍欠采样等距 mask：

%% Displaying mask.
% This code functions to check whether the mask has been made correctly.
% White lines indicate 1's. Black lines indicate 0's. All values should be either 0 or 1.
% There should be a white band in the middle, with striped lines surrounding it.

% 取均值其实没什么特殊意义, 只是把 5 线圈对应的 5 通道压缩成单通道而已
% dispMask = mask(:, :, 1);  % 用 mask(:,:,1) 单独取出一个通道也行, 2,3,4,5 均亦可

dispMask = mean(mask, 3);  % M = mean(A,dim) returns the mean along dimension dim
figure;
imagesc(dispMask);
colormap('gray');
colorbar();
title('equispaced mask (R = 5)')

7. 模拟/生成欠采样 k 空间数据并显示

fullKpace 为各线圈全采样 MR 图像对应的 k 空间域数据，size = 120×128×5，type = complex double

mask 为 0/1 等距掩模，用于模拟/生成 R = 5 倍欠采样 k 空间数据，size = 120×128×5，type = complex double

downSampledKpace 为 fullKpace 和 mask 通过阵列乘法 —— 按元素相乘得到的模拟 R = 5 倍欠采样 k 空间复数值数据，size = 120×128×5，type = complex double

%% Generating the down-sampled k-space image
% Obtain the Hadamard product (elementwise matrix multiplication) between the 
% full k-space data and the mask.
% 
% Hint: use the .* operator, not the * operator, which does matrix multiplication.

% Elementwise (Hadamard) product (multiplication) of matrices.  % 阵列乘法 —— 按元素相乘
% 5 个 fullKspace 各不相同 (位置及敏感度差异导致的观测区别)
downSampledKspace = fullKspace .* mask;  % 用 masked data 模拟欠采样 k-space (复数值) 128×120×5 complex double

downSampleKspaceImage 是通过 平方和根 (root-sum-of-square) 方法 组合所有线圈各自的欠采样 masked k 空间 复数值数据 重建的总体欠采样 masked k 空间 实数值频谱，size = 120×128，type = double

downSampleKspaceImage_Nc 则分别为由各线圈各自的欠采样 masked k 空间 复数值数据 取模得到的 实数值频谱，size = 120×128，type = double

%% Confirmation
% Confirming that the masking has been performed properly in k-space.

% 实数值频谱, 仅用于展示, 无实际用途
downSampledKspaceImage = rsos(downSampledKspace);  % RSOS 组合多个欠采样 k 空间, 重建得到 实值欠采样 k 空间频谱图像 (多对一) 128×120 double
% figure;
% imagesc(downSampledKspaceImage, [0, 10000]);
% colormap('gray');
% colorbar();
assert(isequal(size(downSampledKspace), size(brainCoils)), 'Reconstruction is of the wrong shape.')

figure;
colormap('gray');

subplot(2, 3, 1);
downSampledKspaceImage_1 = rsos(downSampledKspace(:,:,1));
imagesc(downSampledKspaceImage_1, [0, 10000]);
title('coil-1 masked k-space');
colorbar();

subplot(2, 3, 2);
downSampledKspaceImage_2 = rsos(downSampledKspace(:,:,2));
imagesc(downSampledKspaceImage_2, [0, 10000]);
title('coil-2 masked k-space');
colorbar();

subplot(2, 3, 3);
downSampledKspaceImage_3 = rsos(downSampledKspace(:,:,3));
imagesc(downSampledKspaceImage_3, [0, 10000]);
title('coil-3 masked k-space');
colorbar();

subplot(2, 3, 4);
downSampledKspaceImage_4 = rsos(downSampledKspace(:,:,4));
imagesc(downSampledKspaceImage_4, [0, 10000]);
title('coil-4 masked k-space');
colorbar();

subplot(2, 3, 5);
downSampledKspaceImage_5 = rsos(downSampledKspace(:,:,5));
imagesc(downSampledKspaceImage_5, [0, 10000]);
title('coil-5 masked k-space');
colorbar();

subplot(2, 3, 6);
imagesc(downSampledKspaceImage, [0, 10000]);
title('RSOS reconstruction masked k-space');
colorbar();

8. 欠采样 k 空间数据 IFFT 转换到图像域并显示

dsBrainCoils 为 R = 5 倍欠采样 k 空间 复数值数据 的 downSampleKspace 图像域 MR 图像 (Nc = 5 个线圈各自的都叠在一起)，size = 120×128×5，type = complex double

%% Returning the downsampled data to image space.
%% complex sub-sampled MR image

dsBrainCoils = ifftshift(ifft2(fftshift(downSampledKspace)));  % 注意转换到图像域时, IFFT 前后要分别 shift 120x128x5 complex double
assert(isequal(size(dsBrainCoils), size(brainCoils)), 'Image shape is wrong.');

dsImage 是通过 平方和根 (root-sum-of-square) 方法 组合所有线圈各自的欠采样 MR 复数值图像数据 dsBrainCoils 而重建的总体欠采样 MR 实数值图像，size = 120×128，type = double

dsImage_Nc 则分别为由各线圈各自的欠采样 MR 复数值图像数据取模得到的 实数值 MR 图像，size = 120×128，type = double

%% Assignment 2-1: SENSE Reconstruction
%% SENSE
% (SENSitivity Encoding) parallel imaging reconstruction method.
% 
% See <http://mriquestions.com/senseasset.html here> for a basic introduction. 
% See <https://www.researchgate.net/publication/8178233_SMASH_SENSE_PILS_GRAPPA_how_to_choose_the_optimal_method 
% this review paper> (p.226) for a more detailed explanation.
% 
% The notation used in the code below is based on the review paper.

%% Loading data
brainCoilsData = load('brain_coil.mat');
brainCoils = brainCoilsData.brain_coil_tmp;  % 120×128×5 double

coilMapData = load('coil_sensitivity_map.mat');
coilMap = coilMapData.coil_map;  % 120×128×5 double

%% Examine full sampled MR images of each coil
figure;
colormap('gray');

subplot(2, 3, 1);
imagesc(brainCoils(:,:,1));
title('coil-1 full sampled MR image');
colorbar();

subplot(2, 3, 2);
imagesc(brainCoils(:,:,2));
title('coil-2 full sampled MR image');
colorbar();

subplot(2, 3, 3);
imagesc(brainCoils(:,:,3));
title('coil-3 full sampled MR image');
colorbar();

subplot(2, 3, 4);
imagesc(brainCoils(:,:,4));
title('coil-4 full sampled MR image');
colorbar();

subplot(2, 3, 5);
imagesc(brainCoils(:,:,5));
title('coil-5 full sampled MR image');
colorbar();

%% Examine sensitivity map of each coil
figure;
colormap('gray');

subplot(2, 3, 1);
imagesc(coilMap(:,:,1));
title('coil-1 sensitivity map');
colorbar();

subplot(2, 3, 2);
imagesc(coilMap(:,:,2));
title('coil-2 sensitivity map');
colorbar();

subplot(2, 3, 3);
imagesc(coilMap(:,:,3));
title('coil-3 sensitivity map');
colorbar();

subplot(2, 3, 4);
imagesc(coilMap(:,:,4));
title('coil-4 sensitivity map');
colorbar();

subplot(2, 3, 5);
imagesc(coilMap(:,:,5));
title('coil-5 sensitivity map');
colorbar();

%% test about generating sensitivity map manually
BodybrainCoils = rsos(brainCoils);  % RSOS reconstruction for body MR image

figure;
colormap('gray');

subplot(2, 3, 1);
sensitivity_map_1 = brainCoils(:,:,1) ./ BodybrainCoils;
imagesc(sensitivity_map_1);
title('coil-1 sensitivity map (manual)');
colorbar();

subplot(2, 3, 2);
sensitivity_map_2 = brainCoils(:,:,2) ./ BodybrainCoils;
imagesc(sensitivity_map_2);
title('coil-2 sensitivity map (manual)');
colorbar();

subplot(2, 3, 3);
sensitivity_map_3 = brainCoils(:,:,3) ./ BodybrainCoils;
imagesc(sensitivity_map_3);
title('coil-3 sensitivity map (manual)');
colorbar();

subplot(2, 3, 4);
sensitivity_map_4 = brainCoils(:,:,4) ./ BodybrainCoils;
imagesc(sensitivity_map_4);
title('coil-4 sensitivity map (manual)');
colorbar();

subplot(2, 3, 5);
sensitivity_map_5 = brainCoils(:,:,5) ./ BodybrainCoils;
imagesc(sensitivity_map_5);
title('coil-5 sensitivity map (manual)');
colorbar();

subplot(2, 3, 6);
imagesc(BodybrainCoils);
title('RSOS reconstruction body MR image');
colorbar();

%% size
[phaseLength, freqLength, coilNum] = size(brainCoils);  % 120, 128, Nc=5
disp(size(brainCoils));  % Phases, Frequencies, Coil Number

%% Convert to Fourier Domain with 2D Fast Fourier Transform
% Remember to use fftshift and ifftshift when using the Fourier Transform. 
% Otherwise, the k-space will not be centered properly.

fullKspace = ifftshift(fft2(fftshift(brainCoils)));  % 注意转换到 k 空间时, FFT 前后要分别 shift 120×128×5 complex double

%% Examine raw k-space.
fullKspaceImage = rsos(fullKspace);  % RSOS 重建得到实值全采样 k 空间频谱图像

figure;
colormap('gray');

subplot(2, 3, 1);
imagesc(rsos(fullKspace(:,:,1)), [0, 10000]);
title('coil-1 k-space');
colorbar();

subplot(2, 3, 2);
imagesc(rsos(fullKspace(:,:,2)), [0, 10000]);
title('coil-2 k-space');
colorbar();

subplot(2, 3, 3);
imagesc(rsos(fullKspace(:,:,3)), [0, 10000]);
title('coil-3 k-space');
colorbar();

subplot(2, 3, 4);
imagesc(rsos(fullKspace(:,:,4)), [0, 10000]);
title('coil-4 k-space');
colorbar();

subplot(2, 3, 5);
imagesc(rsos(fullKspace(:,:,5)), [0, 10000]);
title('coil-5 k-space');
colorbar();

subplot(2, 3, 6);
imagesc(fullKspaceImage, [0, 10000]);
title('RSOS reconstruction k-space');
colorbar();

%% Setting parameters for later use.
%% setting reduction/accelation factor
downSamplingRate = 5;  % 缩减/加速因子 R = 5 —— 5 倍欠采样/加速

%% Making the mask
% Creating a mask for the downsampling and the ACS lines.

mask = zeros(size(fullKspace));  % Making a mask for the brain coils.
for line = 1:phaseLength   % Goes along the phase encoding axis. (vertical direction ky - dim=0)
    if rem(line, downSamplingRate) == 0  % r = rem(a,b) returns the remainder after division of a by b, 
        mask(line, :, :) = 1;  % Broadcasting the value of 1 to mask(i, :, :).  % 每 R 行采集一条相位编码线, mask=1
    end    
end
assert(isequal(size(mask), size(brainCoils)), 'Mask size is incorrect.')  % mask 要和 full FOV MR image 等尺寸

%% Displaying mask.
% This code functions to check whether the mask has been made correctly.
% White lines indicate 1's. Black lines indicate 0's. All values should be either 0 or 1.
% There should be a white band in the middle, with striped lines surrounding it.

% 取均值其实没什么特殊意义, 只是把 5 线圈对应的 5 通道压缩成单通道而已
% dispMask = mask(:, :, 1);  % 用 mask(:,:,1) 单独取出一个通道也行, 2,3,4,5 均亦可

dispMask = mean(mask, 3);  % M = mean(A,dim) returns the mean along dimension dim
figure;
imagesc(dispMask);
colormap('gray');
colorbar();
title('equispaced mask (R = 5)')

%% Generating the down-sampled k-space image
% Obtain the Hadamard product (elementwise matrix multiplication) between the 
% full k-space data and the mask.
% 
% Hint: use the .* operator, not the * operator, which does matrix multiplication.

% Elementwise (Hadamard) product (multiplication) of matrices.  % 阵列乘法 —— 按元素相乘
% 5 个 fullKspace 各不相同 (位置及敏感度差异导致的观测区别)
downSampledKspace = fullKspace .* mask;  % 用 masked data 模拟欠采样 k-space (复数值) 128×120×5 complex double

%% Confirmation
% Confirming that the masking has been performed properly in k-space.

% 实数值频谱, 仅用于展示, 无实际用途
downSampledKspaceImage = rsos(downSampledKspace);  % RSOS 组合多个欠采样 k 空间, 重建得到 实值欠采样 k 空间频谱图像 (多对一) 128×120 double
% figure;
% imagesc(downSampledKspaceImage, [0, 10000]);
% colormap('gray');
% colorbar();
assert(isequal(size(downSampledKspace), size(brainCoils)), 'Reconstruction is of the wrong shape.')

figure;
colormap('gray');

subplot(2, 3, 1);
downSampledKspaceImage_1 = rsos(downSampledKspace(:,:,1));
imagesc(downSampledKspaceImage_1, [0, 10000]);
title('coil-1 masked k-space');
colorbar();

subplot(2, 3, 2);
downSampledKspaceImage_2 = rsos(downSampledKspace(:,:,2));
imagesc(downSampledKspaceImage_2, [0, 10000]);
title('coil-2 masked k-space');
colorbar();

subplot(2, 3, 3);
downSampledKspaceImage_3 = rsos(downSampledKspace(:,:,3));
imagesc(downSampledKspaceImage_3, [0, 10000]);
title('coil-3 masked k-space');
colorbar();

subplot(2, 3, 4);
downSampledKspaceImage_4 = rsos(downSampledKspace(:,:,4));
imagesc(downSampledKspaceImage_4, [0, 10000]);
title('coil-4 masked k-space');
colorbar();

subplot(2, 3, 5);
downSampledKspaceImage_5 = rsos(downSampledKspace(:,:,5));
imagesc(downSampledKspaceImage_5, [0, 10000]);
title('coil-5 masked k-space');
colorbar();

subplot(2, 3, 6);
imagesc(downSampledKspaceImage, [0, 10000]);
title('RSOS reconstruction masked k-space');
colorbar();


%% Returning the downsampled data to image space.
%% complex sub-sampled MR image
dsBrainCoils = ifftshift(ifft2(fftshift(downSampledKspace)));  % 注意转换到图像域时, IFFT 前后要分别 shift 120x128x5 complex double
assert(isequal(size(dsBrainCoils), size(brainCoils)), 'Image shape is wrong.');

%% Checking reconstructed image.
%%
% Using rsos (root sum of squares) to make visualization possible.
dsImage = rsos(dsBrainCoils);  % RSOS 重建得到 实值欠采样空间域混叠 MR 图像 120×128 double
% figure;
% imagesc(dsImage);
% colormap('gray');
% colorbar();

figure;
colormap('gray');

subplot(2, 3, 1);
dsImage_1 = rsos(dsBrainCoils(:,:,1)) ;
imagesc(dsImage_1);
colormap('gray');
colorbar();
title('coil-1 partial FOV MR image');

subplot(2, 3, 2);
dsImage_2 = rsos(dsBrainCoils(:,:,2)) ;
imagesc(dsImage_2);
colormap('gray');
colorbar();
title('coil-2 partial FOV MR image');

subplot(2, 3, 3);
dsImage_3 = rsos(dsBrainCoils(:,:,3)) ;
imagesc(dsImage_3);
colormap('gray');
colorbar();
title('coil-3 partial FOV MR image');

subplot(2, 3, 4);
dsImage_4 = rsos(dsBrainCoils(:,:,4)) ;
imagesc(dsImage_4);
colormap('gray');
colorbar();
title('coil-4 partial FOV MR image');

subplot(2, 3, 5);
dsImage_5 = rsos(dsBrainCoils(:,:,5)) ;
imagesc(dsImage_5);
colormap('gray');
colorbar();
title('coil-5 partial FOV MR image');

subplot(2, 3, 6);
imagesc(dsImage);
title('RSOS reconstruction partial FOV MR image');
colorbar();


%% SENSE Reconstruction
%%
% Not a perfect solution but works for this experiment.
fieldOfView = floor(phaseLength/downSamplingRate);  % FOV_y = 120 / 5 = 24 
                                                    % FOV_x = 128
senseRecon = zeros(phaseLength, freqLength);  % 120×128 complex double - MR 重建结果
vectorI = zeros(coilNum, 1);  % 5×1 complex double - 混叠像素向量(5个线圈在同一像素位置的复值)
cHat = zeros(coilNum, downSamplingRate);  % 5×5 double - 敏感度图(线圈数×重叠像素数)

for x = 1:freqLength  % kx = FOV_x = 128
    for y = 1:fieldOfView  % FOV_y = 24
        % 欠采样 MR 图像在图像域 pixel(x, y) 处的 Nc=5 个混叠像素值存入 5×1 complex double 向量 I
        vectorI = reshape(dsBrainCoils(y, x, :), coilNum, 1);  % All channels of single pixel in image.
        % 每个混叠像素 pixel(x, y) 都由 R = 5 个 desired full MR image 的像素重叠而成
        for k = 1:downSamplingRate  % R = 5
            % Coil sensitivities of all channels of a single pixel per segment.
            cHat(:, k) = reshape(coilMap(y + (k-1) .* fieldOfView, x, :), coilNum, 1);  % 5×1 double
                % y + (k-1).*fieldOfView 是敏感度图 FOV_y 方向上的偏移量 
                % cHat(:, k) 表示所有 5 个线圈在第 k 个 desired full MR image 像素点处的敏感度值
        end
        
        % Eq.4
        cHatPinv = pinv(cHat);  % Moore-Penrose Pseudoinverse.
        vectorRho = cHatPinv * vectorI;  % 5×1 complex double - 5 个 desired full MR image 像素点
        
        for k = 1:downSamplingRate
            % 按 FOV_y 方向上的偏移量等距保存 desired full MR image 像素点
            senseRecon(y + (k-1) .* fieldOfView, x) = vectorRho(k);
        end
    end
end

%% Visualizing results
%%
senseImage = rsos(senseRecon);
senseImage = senseImage .* downSamplingRate;  % Necessary to correct for the missing values.
original = rsos(brainCoils);
delta = original-senseImage;

figure;
colormap('gray');

subplot(2, 2, 1);
imagesc(original);
title('Original Image');
colorbar();

subplot(2, 2, 2);
imagesc(dsImage);
title('Downsampled Image');
colorbar();

subplot(2, 2, 3);
imagesc(senseImage);
title('SENSE Image');
colorbar();

subplot(2, 2, 4);
imagesc(delta);
title('Difference Image');
colorbar();

为什么欠采样 k 空间在 ky 方向间隔增大 5 倍，IFFT 到图像域时却非 1/5 FOV_y 的 MR 图像？

注意，根据《数字信号处理》理论，频域中的有限离散值，在时域中对应的有限长序列为：原序列以采样点数为周期进行周期延拓后的主值序列。

因此，此处每张混叠部分 FOV MR 图像显示的是长度为 120 / 5 = 24 的主值序列的 R = 5 次周期延拓的结果，因此呈现的混叠部分 FOV MR 图像视野高度 FOV_y 仍为 120。若要得到理论上的长度为 24 的混叠部分 FOV MR 图像，只需任取 R = 5 个周期中的任一者即可。例如，对于 RSOS 重建的混叠部分 FOV MR 图像 (最后一张)，其理论形象为 (1/5 FOV_y)：

9. SENSE 重建

9.1 原理

9.2 实现

%% SENSE Reconstruction
%%
% Not a perfect solution but works for this experiment.
fieldOfView = floor(phaseLength/downSamplingRate);  % FOV_y = 120 / 5 = 24 
                                                    % FOV_x = 128
senseRecon = zeros(phaseLength, freqLength);  % 120×128 complex double - MR 重建结果
vectorI = zeros(coilNum, 1);  % 5×1 complex double - 混叠像素向量(5个线圈在同一像素位置的复值)
cHat = zeros(coilNum, downSamplingRate);  % 5×5 double - 敏感度图(线圈数×重叠像素数)

for x = 1:freqLength  % kx = FOV_x = 128
    for y = 1:fieldOfView  % FOV_y = 24
        % 欠采样 MR 图像在图像域 pixel(x, y) 处的 Nc=5 个混叠像素值存入 5×1 complex double 向量 I
        vectorI = reshape(dsBrainCoils(y, x, :), coilNum, 1);  % All channels of single pixel in image.
        % 每个混叠像素 pixel(x, y) 都由 R = 5 个 desired full MR image 的像素重叠而成
        for k = 1:downSamplingRate  % R = 5
            % Coil sensitivities of all channels of a single pixel per segment.
            cHat(:, k) = reshape(coilMap(y + (k-1) .* fieldOfView, x, :), coilNum, 1);  % 5×1 double
                % y + (k-1).*fieldOfView 是敏感度图 FOV_y 方向上的偏移量 
                % cHat(:, k) 表示所有 5 个线圈在第 k 个 desired full MR image 像素点处的敏感度值
        end
        
        % Eq.4
        cHatPinv = pinv(cHat);  % Moore-Penrose Pseudoinverse. 伪逆
        vectorRho = cHatPinv * vectorI;  % 5×1 complex double - 5 个 desired full MR image 像素点
        
        for k = 1:downSamplingRate
            % 按 FOV_y 方向上的偏移量等距保存 desired full MR image 像素点
            senseRecon(y + (k-1) .* fieldOfView, x) = vectorRho(k);
        end
    end
end

Moore-Penrose 伪逆

Moore-Penrose 伪逆 是一种矩阵，可在不存在逆矩阵的情况下作为逆矩阵的部分替代。此矩阵常被用于求解没有唯一解或有许多解的线性方程组。

对于任何矩阵 A 来说，伪逆 B 都存在，是唯一的，并且具有与 A' (A 的转置) 相同的维度。若 A 是方阵且非奇异，则 pinv(A) 只是一种成本比较高的计算 inv(A) 的方式。但若 A 不是方阵，或是方阵且奇异，则 inv(A) 不存在。在这些情况下，pinv(A) 拥有 inv(A) 的部分（但非全部）属性：

伪逆计算基于 svd(A)。该计算将小于 tol 的奇异值视为零。

具体地：

pinv 通过奇异值分解来形成 A 的伪逆。S 对角线上小于奇异值容差 tol 的奇异值 (不那么重要或代表性不够的特征值) 被视为零，而 A 的表示由此变成：

因此 A 的伪逆等于：

参考文献：Moore-Penrose 伪逆 - MATLAB pinv- MathWorks 中国

9.3 可视化

%% Visualizing results
%%
senseImage = rsos(senseRecon);
senseImage = senseImage .* downSamplingRate;  % Necessary to correct for the missing values.
original = rsos(brainCoils);
delta = original-senseImage;

figure;
colormap('gray');

subplot(2, 2, 1);
imagesc(original);
title('Original Image');
colorbar();

subplot(2, 2, 2);
imagesc(dsImage);
title('Downsampled Image');
colorbar();

subplot(2, 2, 3);
imagesc(senseImage);
title('SENSE Image');
colorbar();

subplot(2, 2, 4);
imagesc(delta);
title('Difference Image');
colorbar();

完整程序 SENSE.m

%% Assignment 2-1: SENSE Reconstruction
%% SENSE
% (SENSitivity Encoding) parallel imaging reconstruction method.
% 
% See <http://mriquestions.com/senseasset.html here> for a basic introduction. 
% See <https://www.researchgate.net/publication/8178233_SMASH_SENSE_PILS_GRAPPA_how_to_choose_the_optimal_method 
% this review paper> (p.226) for a more detailed explanation.
% 
% The notation used in the code below is based on the review paper.

%% Loading data
brainCoilsData = load('brain_coil.mat');
brainCoils = brainCoilsData.brain_coil_tmp;  % 120×128×5 double

coilMapData = load('coil_sensitivity_map.mat');
coilMap = coilMapData.coil_map;  % 120×128×5 double

%% Examine full sampled MR images of each coil
figure;
colormap('gray');

subplot(2, 3, 1);
imagesc(brainCoils(:,:,1));
title('coil-1 full sampled MR image');
colorbar();

subplot(2, 3, 2);
imagesc(brainCoils(:,:,2));
title('coil-2 full sampled MR image');
colorbar();

subplot(2, 3, 3);
imagesc(brainCoils(:,:,3));
title('coil-3 full sampled MR image');
colorbar();

subplot(2, 3, 4);
imagesc(brainCoils(:,:,4));
title('coil-4 full sampled MR image');
colorbar();

subplot(2, 3, 5);
imagesc(brainCoils(:,:,5));
title('coil-5 full sampled MR image');
colorbar();

%% Examine sensitivity map of each coil
figure;
colormap('gray');

subplot(2, 3, 1);
imagesc(coilMap(:,:,1));
title('coil-1 sensitivity map');
colorbar();

subplot(2, 3, 2);
imagesc(coilMap(:,:,2));
title('coil-2 sensitivity map');
colorbar();

subplot(2, 3, 3);
imagesc(coilMap(:,:,3));
title('coil-3 sensitivity map');
colorbar();

subplot(2, 3, 4);
imagesc(coilMap(:,:,4));
title('coil-4 sensitivity map');
colorbar();

subplot(2, 3, 5);
imagesc(coilMap(:,:,5));
title('coil-5 sensitivity map');
colorbar();

%% size
[phaseLength, freqLength, coilNum] = size(brainCoils);  % 120, 128, Nc=5
disp(size(brainCoils));  % Phases, Frequencies, Coil Number

%% Convert to Fourier Domain with 2D Fast Fourier Transform
% Remember to use fftshift and ifftshift when using the Fourier Transform. 
% Otherwise, the k-space will not be centered properly.

fullKspace = ifftshift(fft2(fftshift(brainCoils)));  % 注意转换到 k 空间时, FFT 前后要分别 shift 120×128×5 complex double

%% Examine raw k-space.
fullKspaceImage = rsos(fullKspace);  % RSOS 重建得到实值全采样 k 空间频谱图像

figure;
colormap('gray');

subplot(2, 3, 1);
imagesc(rsos(fullKspace(:,:,1)), [0, 10000]);
title('coil-1 k-space');
colorbar();

subplot(2, 3, 2);
imagesc(rsos(fullKspace(:,:,2)), [0, 10000]);
title('coil-2 k-space');
colorbar();

subplot(2, 3, 3);
imagesc(rsos(fullKspace(:,:,3)), [0, 10000]);
title('coil-3 k-space');
colorbar();

subplot(2, 3, 4);
imagesc(rsos(fullKspace(:,:,4)), [0, 10000]);
title('coil-4 k-space');
colorbar();

subplot(2, 3, 5);
imagesc(rsos(fullKspace(:,:,5)), [0, 10000]);
title('coil-5 k-space');
colorbar();

subplot(2, 3, 6);
imagesc(fullKspaceImage, [0, 10000]);
title('RSOS reconstruction k-space');
colorbar();

%% Setting parameters for later use.
%% setting reduction/accelation factor
downSamplingRate = 5;  % 缩减/加速因子 R = 5 —— 5 倍欠采样/加速

%% Making the mask
% Creating a mask for the downsampling and the ACS lines.

mask = zeros(size(fullKspace));  % Making a mask for the brain coils.
for line = 1:phaseLength   % Goes along the phase encoding axis. (vertical direction ky - dim=0)
    if rem(line, downSamplingRate) == 0  % r = rem(a,b) returns the remainder after division of a by b, 
        mask(line, :, :) = 1;  % Broadcasting the value of 1 to mask(i, :, :).  % 每 R 行采集一条相位编码线, mask=1
    end    
end
assert(isequal(size(mask), size(brainCoils)), 'Mask size is incorrect.')  % mask 要和 full FOV MR image 等尺寸

%% Displaying mask.
% This code functions to check whether the mask has been made correctly.
% White lines indicate 1's. Black lines indicate 0's. All values should be either 0 or 1.
% There should be a white band in the middle, with striped lines surrounding it.

% 取均值其实没什么特殊意义, 只是把 5 线圈对应的 5 通道压缩成单通道而已
% dispMask = mask(:, :, 1);  % 用 mask(:,:,1) 单独取出一个通道也行, 2,3,4,5 均亦可

dispMask = mean(mask, 3);  % M = mean(A,dim) returns the mean along dimension dim
figure;
imagesc(dispMask);
colormap('gray');
colorbar();
title('equispaced mask (R = 5)')

%% Generating the down-sampled k-space image
% Obtain the Hadamard product (elementwise matrix multiplication) between the 
% full k-space data and the mask.
% 
% Hint: use the .* operator, not the * operator, which does matrix multiplication.

% Elementwise (Hadamard) product (multiplication) of matrices.  % 阵列乘法 —— 按元素相乘
% 5 个 fullKspace 各不相同 (位置及敏感度差异导致的观测区别)
downSampledKspace = fullKspace .* mask;  % 用 masked data 模拟欠采样 k-space (复数值) 128×120×5 complex double

%% Confirmation
% Confirming that the masking has been performed properly in k-space.

% 实数值频谱, 仅用于展示, 无实际用途
downSampledKspaceImage = rsos(downSampledKspace);  % RSOS 组合多个欠采样 k 空间, 重建得到 实值欠采样 k 空间频谱图像 (多对一) 128×120 double
% figure;
% imagesc(downSampledKspaceImage, [0, 10000]);
% colormap('gray');
% colorbar();
assert(isequal(size(downSampledKspace), size(brainCoils)), 'Reconstruction is of the wrong shape.')

figure;
colormap('gray');

subplot(2, 3, 1);
downSampledKspaceImage_1 = rsos(downSampledKspace(:,:,1));
imagesc(downSampledKspaceImage_1, [0, 10000]);
title('coil-1 masked k-space');
colorbar();

subplot(2, 3, 2);
downSampledKspaceImage_2 = rsos(downSampledKspace(:,:,2));
imagesc(downSampledKspaceImage_2, [0, 10000]);
title('coil-2 masked k-space');
colorbar();

subplot(2, 3, 3);
downSampledKspaceImage_3 = rsos(downSampledKspace(:,:,3));
imagesc(downSampledKspaceImage_3, [0, 10000]);
title('coil-3 masked k-space');
colorbar();

subplot(2, 3, 4);
downSampledKspaceImage_4 = rsos(downSampledKspace(:,:,4));
imagesc(downSampledKspaceImage_4, [0, 10000]);
title('coil-4 masked k-space');
colorbar();

subplot(2, 3, 5);
downSampledKspaceImage_5 = rsos(downSampledKspace(:,:,5));
imagesc(downSampledKspaceImage_5, [0, 10000]);
title('coil-5 masked k-space');
colorbar();

subplot(2, 3, 6);
imagesc(downSampledKspaceImage, [0, 10000]);
title('RSOS reconstruction masked k-space');
colorbar();


%% Returning the downsampled data to image space.
%% complex sub-sampled MR image
dsBrainCoils = ifftshift(ifft2(fftshift(downSampledKspace)));  % 注意转换到图像域时, IFFT 前后要分别 shift 120x128x5 complex double
assert(isequal(size(dsBrainCoils), size(brainCoils)), 'Image shape is wrong.');

%% Checking reconstructed image.
%%
% Using rsos (root sum of squares) to make visualization possible.
dsImage = rsos(dsBrainCoils);  % RSOS 重建得到 实值欠采样空间域混叠 MR 图像 120×128 double
% figure;
% imagesc(dsImage);
% colormap('gray');
% colorbar();

figure;
colormap('gray');

subplot(2, 3, 1);
dsImage_1 = rsos(dsBrainCoils(:,:,1)) ;
imagesc(dsImage_1);
colormap('gray');
colorbar();
title('coil-1 partial FOV MR image');

subplot(2, 3, 2);
dsImage_2 = rsos(dsBrainCoils(:,:,2)) ;
imagesc(dsImage_2);
colormap('gray');
colorbar();
title('coil-2 partial FOV MR image');

subplot(2, 3, 3);
dsImage_3 = rsos(dsBrainCoils(:,:,3)) ;
imagesc(dsImage_3);
colormap('gray');
colorbar();
title('coil-3 partial FOV MR image');

subplot(2, 3, 4);
dsImage_4 = rsos(dsBrainCoils(:,:,4)) ;
imagesc(dsImage_4);
colormap('gray');
colorbar();
title('coil-4 partial FOV MR image');

subplot(2, 3, 5);
dsImage_5 = rsos(dsBrainCoils(:,:,5)) ;
imagesc(dsImage_5);
colormap('gray');
colorbar();
title('coil-5 partial FOV MR image');

subplot(2, 3, 6);
imagesc(dsImage);
title('RSOS reconstruction partial FOV MR image');
colorbar();


%% SENSE Reconstruction
%%
% Not a perfect solution but works for this experiment.
fieldOfView = floor(phaseLength/downSamplingRate);  % FOV_y = 120 / 5 = 24 
                                                    % FOV_x = 128
senseRecon = zeros(phaseLength, freqLength);  % 120×128 complex double - MR 重建结果
vectorI = zeros(coilNum, 1);  % 5×1 complex double - 混叠像素向量(5个线圈在同一像素位置的复值)
cHat = zeros(coilNum, downSamplingRate);  % 5×5 double - 敏感度图(线圈数×重叠像素数)

for x = 1:freqLength  % kx = FOV_x = 128
    for y = 1:fieldOfView  % FOV_y = 24
        % 欠采样 MR 图像在图像域 pixel(x, y) 处的 Nc=5 个混叠像素值存入 5×1 complex double 向量 I
        vectorI = reshape(dsBrainCoils(y, x, :), coilNum, 1);  % All channels of single pixel in image.
        % 每个混叠像素 pixel(x, y) 都由 R = 5 个 desired full MR image 的像素重叠而成
        for k = 1:downSamplingRate  % R = 5
            % Coil sensitivities of all channels of a single pixel per segment.
            cHat(:, k) = reshape(coilMap(y + (k-1) .* fieldOfView, x, :), coilNum, 1);  % 5×1 double
                % y + (k-1).*fieldOfView 是敏感度图 FOV_y 方向上的偏移量 
                % cHat(:, k) 表示所有 5 个线圈在第 k 个 desired full MR image 像素点处的敏感度值
        end
        
        % Eq.4
        cHatPinv = pinv(cHat);  % Moore-Penrose Pseudoinverse.
        vectorRho = cHatPinv * vectorI;  % 5×1 complex double - 5 个 desired full MR image 像素点
        
        for k = 1:downSamplingRate
            % 按 FOV_y 方向上的偏移量等距保存 desired full MR image 像素点
            senseRecon(y + (k-1) .* fieldOfView, x) = vectorRho(k);
        end
    end
end

%% Visualizing results
%%
senseImage = rsos(senseRecon);
senseImage = senseImage .* downSamplingRate;  % Necessary to correct for the missing values.
original = rsos(brainCoils);
delta = original-senseImage;

figure;
colormap('gray');

subplot(2, 2, 1);
imagesc(original);
title('Original Image');
colorbar();

subplot(2, 2, 2);
imagesc(dsImage);
title('Downsampled Image');
colorbar();

subplot(2, 2, 3);
imagesc(senseImage);
title('SENSE Image');
colorbar();

subplot(2, 2, 4);
imagesc(delta);
title('Difference Image');
colorbar();

参考链接：https://github.com/veritas9872/Medical-Imaging-Tutorial