【MRI】GRAPPA (GeneRalized Autocalibrating Partially Parallel Acquisitions) 算法仿真实验与原理剖析 (Matlab 实现)

何处闻韶

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

阅读量4.5k

点赞数 13

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

于 2020-11-04 09:56:34 首次发布

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

版权

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

2 篇文章 10 订阅

订阅专栏

1. 加载全采样 MR 图像并显示

2. 全采样 MR 图像转换为全采样 k 空间并显示

3. 设置部分参数

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

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

6. 欠采样 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

初始化

%% GRAPPA MRI Reconstruction
%% Initializing
%%
tic
close;
clear;

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

1. 加载全采样 MR 图像并显示

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

%% Loading data
%%
brainCoilsData = load('brain_coil.mat');
brainCoils = brainCoilsData.brain_coil_tmp;  % 120×128×5 double
[phaseLength, freqLength, numCoils] = size(brainCoils);  % 120, 128, 5
disp(size(brainCoils));  % Phases, Frequencies, Coils

%% Show original Brain Image
%%
originalImage = rsos(brainCoils);  % 120×128 double

figure;
colormap('gray');
imagesc(originalImage);  
title('Original Image');
colorbar();

2. 全采样 MR 图像转换为全采样 k 空间并显示

fullKsapce 是由全采样 MR 图像 brainCoils 经 FFT 得到的全采样 k 空间数据，size = 120×128×5，type = complex double

%% Convert to Fourier Domain with 2D Fast Fourier Transform
% Don't forget to use fftshift and ifftshift.

%fullKspace = fftshift(fftshift(fft2(brainCoils), 1), 2);  % 源代码
fullKspace = ifftshift(fft2(fftshift(brainCoils)));  % 修正  120×128×5 complex double

%% Check k-space
%%
figure;
imagesc(rsos(fullKspace), [0, 10000]);
title('Original k-space');
colormap('gray');
colorbar();

fftshift 与 ifftshift

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

举个例子：

3. 设置部分参数

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

acsLength 即 ACS 区域长度 / ACS 相位编码线数目，此处表示中央区域连续 12 条相位编码线作为 ACS。

kernelHeight 指 GRAPPA kernel 高度，kernelWidth 指 GRAPPA kernel 宽度

remainder 用于确定已采集线 (mask=1) 的位置 (R=2 时用于区分奇数索引/偶数索引)

%% Setting parameters for later use.
% Due to overlap with the downsampling lines, the actual ACS length may be 1 
% or 2 lines longer than specified.

samplingRate = 2;  % reduction factor / accelation factor R = 2
acsLength = 12;  % 12 条 ACS 线
kernelHeight = 4;  % kernel 高度为 4
kernelWidth = 3;  % kernel 宽度为 3
remainder = 0;  % 已采集线为偶数行

assert(mod(kernelHeight, 2) == 0, 'kernel height should be an even number')
assert(mod(kernelWidth, 2) == 1, 'kernel width should be an odd number')
assert(remainder < samplingRate, ...
    'The remainder must be a non-negative integer smaller than the sampling rate.')

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

acsStart 为 ACS 起始行索引，acsFinish 为 ACS 结束行索引，在此范围内的 mask 值恒为 1 表示全采样。

mask 为与 fullKpace 等 size 和 type 的 0/1 等距掩模 (不是随机掩模)，用于模拟/制作 R = 2 倍欠采样 k 空间数据，其中第 55-66 共 12 行为 ACS 区域。

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

mask = zeros(size(fullKspace));  % Making a mask for the brain coils. 120×128×5 double

acsStart = floor((phaseLength - acsLength) ./ 2);  % ACS 起始行索引 (120 - 12) / 2 = 54
acsFinish = acsStart + acsLength;  % ACS 结束行索引 54 + 12 = 66

for idx = 1:phaseLength
    if (acsStart < idx) && (idx <= acsFinish)
        mask(idx, :, :) = 1;  % ACS 区域为 1 值掩模
    end
    
    if mod(idx, samplingRate) == remainder
        mask(idx, :, :) = 1;  % 已采集线设在偶数索引处 (remainder=1 则在奇数索引)
    end
end

assert(isequal(size(mask), size(brainCoils)), 'Mask size is incorrect.')

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

%% Displaying mask.
% This code functions to check whether the mask has been made correctly.
% 
% All values should be either 0 or 1. The presence of any other value indicates 
% an error.
% 
% There should be a white band in the middle, with striped lines surrounding 
% it.

figure;
imagesc(mean(mask, 3));
title('Mask array');
colormap('gray');
colorbar();

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

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

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

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

%% Generating down-sampled k-space
% Obtain the Hadamard product (elementwise matrix multiplication) between the 
% full k-space data and the mask.

downSampledKspace = fullKspace .* mask;  % 120×128×5 complex double
assert(isequal(size(downSampledKspace), size(brainCoils)), 'Reconstruction is of the wrong shape.')

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

% 实数值频谱, 仅用于展示, 无实际用途
downSampledKspaceImage = rsos(downSampledKspace);  % RSOS 组合多个欠采样 k 空间, 重建得到 实值欠采样 k 空间频谱图像 (多对一) 128×120 double

figure;
imagesc(downSampledKspaceImage, [0, 10000]);
title('Downsampled k-space');
colormap('gray');
colorbar();

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

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

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

%% Check Downsampled Image
%%
%dsImage = ifft2(ifftshift(ifftshift(downSampledKspace, 2), 1));
dsBrainCoils = ifftshift(ifft2(fftshift(downSampledKspace)));  % 120×128×5 complex double
dsImage = rsos(dsBrainCoils);  % 120×128 double

然后，显示总体 R = 2 倍欠采样 MR 图像：

figure;
imagesc(dsImage);
title('Downsampled Image');
colormap('gray');
colorbar();

7. 重新修正真正的 ACS 区域

由于使用了等距 mask 模拟 R = 2 倍的欠采样，在原 ACS 区域 (行索引 55-66) 的基础上，必然多出 1 条已采集的相位编码线与之相邻，因此可以重新修正 ACS 区域使之扩大。

由于 mask 是用于模拟欠采样的，而在实际欠采样 k 空间数据中，我们并不知道 mask 的存在。因此，我们需要手动获取等同于 mask 的数据 finder。

finder 中，包括 ACS 区域 (55-66) 内的已采集数据 (偶数行) 的行索引 (相位编码方向) 处值为 1，未采集数据 (奇数行) 的行索引为 0，size = 120×1，type = logical，例如：

acsFinder 则在原定 ACS 区域行索引 (相位编码方向) 处值为 1，其余行索引处值均为 0，size = 120×1，type = double，例如：

acsLines 则包含了最终实际的 ACS 区域行索引，size = 13。之所以为 13 而非 12 条, 是因为多出的 1 条 (第 54 条) 来偶数位置索引的已采样数据。例如：

matrixHeight 表示 kernel / 权重覆盖范围的高度为 7 (每次拟合/估计纵跨 7 条相位编码线，其中 第 1、3、5、7 条 为已采集数据点 (kernelHeight = 4)，第 4 条的中点 为待拟合/待估计数据点 (kernelWidth = 3))。

acsTrueLength 表示最终实际的 ACS 区域长度为 13 (第 54-66 行)。

%% Find True ACS
% Find the length of the true ACS domain, including lines due to downsampling.

%finder = any(downSampledKspace, [2, 3]);  % 原代码
%assert(isequal(finder, mean(mask, [2, 3])), 'The indices are wrong.')  % 原代码
finder = any(any(downSampledKspace, 2), 3);  % finder.size = 120×1 logical
assert(isequal(finder, mean(mean(mask, 2),3)), 'The indices are wrong.')  % mean.size = 120×1 double

acsFinder = zeros(phaseLength, 1);  % 120×1 double
for idx = 2:phaseLength-1
    if finder(idx-1) == 1 && finder(idx) == 1 && finder(idx+1) == 1  
        acsFinder(idx-1:idx+1) = 1;  % ACS 区域置 1
    end
end

% Getting the idices of the ACS lines.
acsLines = find(acsFinder);  % 13×1 double (54-66)
% 之所以为 13 而非 12 条, 是因为多出的 1 条来自偶数位置索引的已采样数据

% Checking whether the parameters fit.
matrixHeight = samplingRate .* (kernelHeight-1) + 1;  % 2 * (4-1) + 1 = 7  表示 kernel/权重覆盖范围的高度为 7
acsTrueLength = length(acsLines);  % 表示真实的 ACS 线长度为 13 

assert(acsTrueLength >= matrixHeight, 'ACS is too short for kernel params')
assert((acsLength <= acsTrueLength) && (acsTrueLength <= acsLength+2), 'Mask is incorrect')

8. 构造 kernel 权重矩阵

8.1 正常过程

根据矩阵乘法公式：Y = X * w → 权重矩阵：w = pinv(X) * Y

对应代码中的变量名与运算：weights = pinv(inMatrix) * outMatrix;

acsPhases 为 7，表示覆盖范围高度为 7 的 GRAPPA kernel，在相位编码方向上沿 ACS 区域 (54-66 共 13 行) 有 7 个不同的起始行。

numKernels 为 897，表示 GRAPPA kernel 在 ACS 区域滑动的总次数为 897，对应着 897 个待拟合的 ACS 点 y。987 = 7 * 128，其中 7 表示 7 个 kernel 范围起始行，128 表示频率编码轴长度为 128 (有 128 个横坐标)。

kernelSize 为 60，表示每个 GRAPPA kernel 的权重个数为 60。60 = 4 * 3 * 5，其中 4 是 kernel 高度，3 是 kernel 宽度，5 是 kernel 通道数，对应 Nc = 5 个线圈，表明每个线圈上的待拟合/待估计点都是由所有 Nc = 5 个线圈的已采集数据点共同参与的。

outSize 为 5，表示 每次拟合/估计 1 个位置的点，Nc = 5 个线圈在此相同位置的点由矩阵运算同时得出，对应 outMatrix 的列数。

%% Building the weight matrix
% Has equation of X * w = Y in mind. (Eq.11)
%     X = inMatrix, Y = outMatrix, w = weights.
% Each kernel adds one row to the X and Y matrices.

acsPhases = acsTrueLength - matrixHeight + 1;          % 13 - 7 + 1 = 7, kernel 纵向移动次数 / 起始行次序
numKernels = acsPhases .* freqLength;                  % 7 * 128 = 897, kernel 在 ACS 上移动的总次数 (训练对数)
kernelSize = kernelHeight .* kernelWidth .* numCoils;  % 4 * 3 * 5 = 60, kernel 包含的权重个数
outSize = numCoils .* (samplingRate - 1);              % 5 * (2 - 1) = 5, outMatrix 的列数

outMatrix 表示 ACS 中待拟合的数据点矩阵 Y (可视为 训练集 labels / targets)，size = 896×5，type = complex double。其中，896 = 7 * 128，表示 ACS 区域中央的 7 行 (第 54, 55, 56, 57, 58, 59, 60)、每行的 128 个点为用于估计 GRAPPA kernel 的待拟合数据点 y，而 Nc = 5 个线圈在同一坐标位置处的 5 个点由矩阵运算同时得出。

inMatrix 对应 ACS 中用于数据点矩阵 X (可视为 训练集 data)，size = 896×60，type = complex double。其中，896 = 7 * 128，表示共有 7 * 128 个待拟合数据点 x；60 = 4 * 3 * 5，表示拟合每个点需要的 ACS 数据点 / 权重个数为 60 个。换言之，每 60 个 ACS 数据点 x / 权重拟合一个 y，共需要 896 组。

% Eq.11 S represents the collected signal in each element coil at some position
inMatrix = zeros(numKernels, kernelSize);  % 896×60 complex double
% Eq.11 S(m) 
outMatrix = zeros(numKernels, outSize);    % 896×5 complex double

hkw 为 GRAPPA kernel 宽度的一半 1，hkh 为 GRAPPA kernel 高度的一半 2。

kidx 作为 kernel 计数器。

hkw = floor(kernelWidth/2);     % Half kernel width.  1
hkh = kernelHeight/2;           % Half kernel height. 2
kidx = 1;  % "Kernel index" for counting the number of kernels.

然后，遍历 7 行 (54-60) 128 列 (1-128) 以拟合共 7 * 128 = 896 个点矩阵 Y。

注意，每次估计 kernel 的覆盖范围是 7 行 3 列，但只有其中 4 行属于 x (如第一次为 54、56、58、60)，列选取则是循环移位的 (如第一次为 128、1、2)，那么待拟合点 y 就在正中心 (如第一次为 (57, 1))。以此类推：

第二次行编号为 54、56、58、60，列编号为 1、2、3，待拟合点 y 在 (57，2)；

...

倒数第二次 (第891次) 行编号为 60、62、64、66，列编号为 126、127、128，待拟合点 y 在 (61，127)；

最后一次 (第892次) 行编号为 60、62、64、66，列编号为 127、128、1，待拟合点 y 在 (61，128)。

% acsline = 54,55,56,57,58,59,60, 表示 kernel 的起始行
for acsLine = acsLines(1:acsPhases)'  % 遍历第 1,2,3,4,5,6,7 位   (acsLines 整体范围从 54 到 66)
    % 第一个 phases = linsapce(54, 54+7-1, 4) = linspace(54, 60, 4) = 54  56  58  60
    phases = linspace(acsLine, acsLine+matrixHeight-1, kernelHeight);  % Phases of the kernel
    % 频率编码轴长度 freqLength = 128, freq 表示当前 kernel 中点的横坐标/频率值 (1-128)
    for freq = 1:freqLength
        % 第一个 freqs = linspace(1-1, 1+1, 3) = linspace(0, 2, 3) = 0  1  2
        freqs = linspace(freq-hkw, freq+hkw, kernelWidth);  % Frequencies of the kernel.
        % 第一个 freqs = mod((0  1  2) - 1, 128) + 1 = 128  1  2
        freqs = mod(freqs-1, freqLength) + 1;  % For circular indexing.
        
        % 第一个 selected = linspace(56+1, 58-1, 2-1) = 57
        selected = linspace(phases(hkh)+1, phases(hkh+1)-1, samplingRate-1);
        % 第一个 selected = mod(57-1, 120) + 1 = 57
        selected = mod(selected-1, phaseLength) + 1;  % Selected Y phases.
        
        % 第一个 tempX = downSampledKspace((54  56  58  60), (128  1  2), :) 
        % 共 4 * 3 * 5 = 60 个数据点
        tempX = downSampledKspace(phases, freqs, :);  % 4×3×5 complex double
        % 第一个 tempY = downSampledKspace(57, 1, :)
        % 共 1 * 1 * 5 = 5 个待拟合点
        tempY = downSampledKspace(selected, freq, :);  % 1×1×5 complex double
        
        % Filling in the matrices row by row.
        inMatrix(kidx, :) = reshape(tempX, 1, kernelSize);  % ACS 训练数据点排在第 kids 行上, size = 1×60
        outMatrix(kidx, :) = reshape(tempY, 1, outSize);    % ACS 待拟合数据点排在第 kids 行上, size = 1×5
        
        kidx = kidx + 1;  % 下一保存行, 共 7 * 128 = 896 行, 对应 896 个 ACS 待拟合数据点
    end
end

assert(isequal(size(inMatrix), [numKernels, kernelSize]), 'Incorrect X matrix');
assert(isequal(size(outMatrix), [numKernels, outSize]), 'Incorrect Y matrix');

最后，使用 Moore-Penrose 伪逆 求解权重矩阵 weights，size = 60×5，type = comples double

% Has equation of X * w = Y in mind. (Eq.11)
%     X = inMatrix, Y = outMatrix, w = weights.
% Each kernel adds one row to the X and Y matrices.

% Eq.11 求解 n^(m) (weights) 
%     pinv(inMatrix)    60×896    complex double
%     outMatrix         896×5     complex double
weights = pinv(inMatrix) * outMatrix;  % Calculate the weight matrix.  60×5 complex double 
%     B = pinv(A) returns the Moore-Penrose pseudoinverse of A (伪逆)
                       
assert(isequal(size(weights), [kernelSize, outSize]), 'Incorrect weight size');

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 中国

8.2 可视化

%% 可视化三个矩阵
figure;
rsos_inMatrix = rsos(inMatrix);
imagesc(rsos_inMatrix, [0, 10000]);  % 896×60 double
title('modulus of X');
colormap('gray');
colorbar();

figure;
rsos_outMatrix = rsos(outMatrix);
imagesc(rsos_outMatrix, [0, 10000]);  % 896×5 double
title('modulus of Y');
colormap('gray');
colorbar();

figure;
rsos_weights = rsos(weights);
imagesc(rsos_weights);  % 60×5 double
title('modulus of weights');
colormap('gray');
colorbar();

9. GARPPA 重建

9.1 原理

9.2 实现

说明

%% GRAPPA Reconstruction
% Performing a naive reconstruction according to first principles causes an 
% overlap problem.
% 
% The lines immediately before and after the ACS lines are not necessarily 
% spaced with the sampling rate as the spacing.
% 
% This causes alteration of the original data if GRAPPA reconstruction is 
% performed naively.
% 
% The solution is to perform reconstruction on a blank, and then overwrite 
% all vlaues with the original data.
% 
% This alleviates the problem of having to do special operations for the 
% values at the edges.
% 
% Also, the lines of k-space at the start or end of k-space may be neglected 
% (depending on the implementation).
% 
% This requires shifting the finder by the sampling rate to look at the phase 
% from one step above.
% 
% If the downsampling does not match due to incorrect k-space dimensions, 
% errors will be overwritten by the final correction process.

% 根据第一个原理进行朴素的重建会导致重叠问题. 
% ACS 线前后的线不一定要以采样率作为间隔来间隔。 
% 如果朴素地执行 GRAPPA 重建，将导致原始数据的更改。 
% 解决方案是对空白重建，然后用原始数据覆盖所有值。    (数据一致性)
% 这减轻了必须对边缘值进行特殊运算的问题。 
% 而且，可以忽略在 k 空间的开始或结尾处的 k 空间的线（取决于实现方式）。 
% 这需要将 finder 按采样率移动，以从上方一步查看相位。 
% 如果由于不正确的 k 空间尺寸导致下采样不匹配，则最终的校正过程将覆写错误。

先确定未采集数据行索引 fillFinder，再逐点通过 GRAPPA kernel 卷积 (矩阵运算形式) 估计缺失点。对于已采集数据点，不经估计直接覆盖使用，保持 数据一致性。换言之，已采集数据点直接用上，未采集数据点才进行估计。

% Find the indices to fill, including at the beginning of k-space.
temp1 = find(diff(finder) == -1) + 1;  % Gets nearly all the lines.
% 待填充行索引???      3、5 ... 51、53、67、69 ... 117、119 共 53 行
%   diff(finder) 在 54-65 行为 0, 其余奇数行为 1, 偶数行为 -1  119×1 double
%   find(diff(finder) == -1) + 1 找到值为 -1 的偶数行的索引, 然后将索引+1,  53×1 double

temp2 = find(diff(circshift(finder, samplingRate)) == -1) - samplingRate + 1;
% 待填充行索引???      1、3 ... 51、53、67、69 ... 115、117 共 53 行
%   circshift(finder, samplingRate) 偶数行及 56-68 行为 1, 其余行为 0, 120×1 logical
%   diff(circshift(finder, samplingRate)) 在 56-67 行为 0, 其余奇数行为 1, 偶数行为 -1, 119×1 double  
%   find(diff(circshift(finder, samplingRate)) == -1) - samplingRate + 1 找到值为 -1 的偶数行的索引, 
%      然后将索引 - 2 + 1 = -1, 53×1 double
                                                                   
% Second line exists to catch the first few empty lines, if any are present.
% 待重建从而填充的缺失行索引 (mask=0)
fillFinder = unique([temp1; temp2]);  % 1、3 ... 51、53、67、69 ... 117、119 共 54 行

% Shift from first fill line to beginning of kernel data.
upShift = (hkh-1) .* samplingRate + 1;  % (2-1) * 2 + 1 = 3
% Shift from first fill line to end of kernel data.
downShift = hkh .* samplingRate - 1;  % 2 * 2 - 1 = 3

% GRAPPA 重建 k 空间
grappaKspace = zeros(size(downSampledKspace));  % 120×128×5 complex double

% 遍历待重建填充的缺失行索引 (mask=0) 1×53 double
for phase = fillFinder'
    % 第一个 phases = linspace(1-3, 1+3, 4) = -2  0  2  4
    phases = linspace(phase-upShift, phase+downShift, kernelHeight);
    % 第一个 phases = mod((-2  0  2  4) - 1, 120) + 1 = 118  120  2  4
    phases = mod(phases-1, phaseLength) + 1;  % Circularly indexed phases.
    % 频率编码轴长度 freqLength = 128, freq 表示当前 kernel 中点的横坐标/频率值 (1-128)
    for freq = 1:freqLength    
        % 第一个 freqs = linspace(1-1, 1+1, 3) = linspace(0, 2, 3) = 0  1  2
        freqs = linspace(freq-hkw, freq+hkw, kernelWidth);  % Frequencies of the kernel.
        % 第一个 freqs = mod((0  1  2) - 1, 128) + 1 = 128  1  2
        freqs = mod(freqs-1, freqLength) + 1;  % Circularly indexed frequencies.   
        
        kernel = downSampledKspace(phases, freqs, :);  % 取出已知数据点, 4×3×5 complex double
        % One line of the input matrix.
        tempX = reshape(kernel, 1, kernelSize);  % 1×60 complex double
        
        % One line of the output matrix.  
        tempY = tempX * weights;   % 加权线性组合得到缺失点, 1×5 complex double
        tempY = reshape(tempY, (samplingRate-1), 1, numCoils);  % 1×1×5 complex double
        
        % Selected lines of the output matrix to be filled in.
        % 第一个 selected = linspace(0+1, 2-1, 2-1) = 1
        selected = linspace(phases(hkh)+1, phases(hkh+1)-1, samplingRate-1);
         % 第一个 selected = mod(1-1, 120) + 1 = 1
        selected = mod(selected-1, phaseLength) + 1;  % Selected Y phases.
        
        grappaKspace(selected, freq, :) = tempY;  % 缺失点估计值填充
    end
end

% Filling in all the original data.
% Doing it this way solves the edge overlap problem.
% 只估计未采集的数据, 而直接拷贝已采集的数据, 从而实现数据一致性 (data consistency)
grappaKspace(finder, :, :) = downSampledKspace(finder, :, :);  % 数据一致性处理

assert(isequal(size(downSampledKspace), size(grappaKspace)), 'Incorrect matrix size.');
% Somewhat redundant assertion now.
% assert(isequal(grappaKspace(finder, :, :), downSampledKspace(finder, :, :)));

9.3 可视化

显示 RSOS 重建 MR 图像：

%% Display recon image
%%
%recon = ifft2(ifftshift(ifftshift(grappaKspace, 2), 1));  源代码
recon = ifftshift(ifft2(fftshift(grappaKspace)));  % 120×128×5 complex double 修正
reconImage = rsos(recon);  % 120×128 double  RSOS 组合重建结果

显示恢复的各线圈及 RSOS 重建的相位角图像：

%% Display phase angle
%temp = recon(:, :, 1);  % 120×128 complex double
% angle(reconImage) 则由于 reconImage 是实数图像, 相位角已经为 0 了

% P = angle(Z) returns the phase angles, in radians, for each element of complex array Z.
% The angles lie between ±π.
% For complex Z, the magnitude R and phase angle theta are given by
%       R = abs(Z)
%       theta = angle(Z)
% and the statement
%       Z = R.*exp(i*theta)
% converts back to the original complex Z.

figure;
colormap('gray');
 
subplot(2, 3, 1);
imagesc(angle(recon(:, :, 1)));
title('coil-1 recon phase angle');
colorbar();
 
subplot(2, 3, 2);
imagesc(angle(recon(:, :, 2)));
title('coil-2 recon phase angle');
colorbar();
 
subplot(2, 3, 3);
imagesc(angle(recon(:, :, 3)));
title('coil-3 recon phase angle');
colorbar();
 
subplot(2, 3, 4);
imagesc(angle(recon(:, :, 4)));
title('coil-4 recon phase angle');
colorbar();

subplot(2, 3, 5);
imagesc(angle(recon(:, :, 5)));
title('coil-5 recon phase angle');
colorbar();
 
subplot(2, 3, 6);
imagesc(angle(reconImage));
title('RSOS recon MR phase angle');
colorbar();

对欠采样重建图像与全采样重建图像作差并显示：

%% Compute difference image
%%
deltaImage = reconImage - originalImage;  % difference image  重建图像与标签图像之差

figure;
imagesc(deltaImage);   
title('Difference Image');
colormap('gray');
axis('off');
colorbar();

汇总对比显示：

%% Summary
%%
figure;

subplot(2, 2, 1);
imagesc(originalImage);
title('Original Image');
colormap('gray');
axis('off');
colorbar();

subplot(2, 2, 2);
imagesc(dsImage);
title('Down-Sampled Image');
colormap('gray');
axis('off');
colorbar();

subplot(2, 2, 3);
imagesc(reconImage);
title('Reconstructed Image');
colormap('gray');
axis('off');
colorbar();

subplot(2, 2, 4);
imagesc(deltaImage);
title('Difference Image');
colormap('gray');
axis('off');
colorbar();

toc

完整程序 GRAPPA.m

%% GRAPPA MRI Reconstruction
%% Initializing
%%
tic
close;
clear;

%% Loading data
%%
brainCoilsData = load('brain_coil.mat');
brainCoils = brainCoilsData.brain_coil_tmp;  % 120×128×5 double
[phaseLength, freqLength, numCoils] = size(brainCoils);  % 120, 128, 5
disp(size(brainCoils));  % Phases, Frequencies, Coils

%% Show original Brain Image
%%
originalImage = rsos(brainCoils);  % 120×128 double

figure;
colormap('gray');
imagesc(originalImage);  
title('Original Image');
colorbar();

%% Convert to Fourier Domain with 2D Fast Fourier Transform
% Don't forget to use fftshift and ifftshift.

%fullKspace = fftshift(fftshift(fft2(brainCoils), 1), 2);  % 源代码
fullKspace = ifftshift(fft2(fftshift(brainCoils)));  % 修正  120×128×5 complex double

%% Check k-space
%%
figure;
imagesc(rsos(fullKspace), [0, 10000]);
title('Original k-space');
colormap('gray');
colorbar();

%% Setting parameters for later use.
% Due to overlap with the downsampling lines, the actual ACS length may be 1 
% or 2 lines longer than specified.

samplingRate = 2;  % reduction factor / accelation factor R = 2
acsLength = 12;  % 12 条 ACS 线
kernelHeight = 4;  % kernel 高度为 4
kernelWidth = 3;  % kernel 宽度为 3
remainder = 0;  % 已采集线为偶数行

assert(mod(kernelHeight, 2) == 0, 'kernel height should be an even number')
assert(mod(kernelWidth, 2) == 1, 'kernel width should be an odd number')
assert(remainder < samplingRate, ...
    'The remainder must be a non-negative integer smaller than the sampling rate.')

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

mask = zeros(size(fullKspace));  % Making a mask for the brain coils. 120×128×5 double

acsStart = floor((phaseLength - acsLength) ./ 2);  % ACS 起始行索引 (120 - 12) / 2 = 54
acsFinish = acsStart + acsLength;  % ACS 结束行索引 54 + 12 = 66

for idx = 1:phaseLength
    if (acsStart < idx) && (idx <= acsFinish)
        mask(idx, :, :) = 1;  % ACS 区域为 1 值掩模
    end
    
    if mod(idx, samplingRate) == remainder
        mask(idx, :, :) = 1;  % 已采集线设在偶数索引处 (remainder=1 则在奇数索引)
    end
end

assert(isequal(size(mask), size(brainCoils)), 'Mask size is incorrect.')

%% Displaying mask.
% This code functions to check whether the mask has been made correctly.
% 
% All values should be either 0 or 1. The presence of any other value indicates 
% an error.
% 
% There should be a white band in the middle, with striped lines surrounding 
% it.

figure;
imagesc(mean(mask, 3));
title('Mask array');
colormap('gray');
colorbar();

%% Generating down-sampled k-space
% Obtain the Hadamard product (elementwise matrix multiplication) between the 
% full k-space data and the mask.

downSampledKspace = fullKspace .* mask;  % 120×128×5 complex double
assert(isequal(size(downSampledKspace), size(brainCoils)), 'Reconstruction is of the wrong shape.')

% 实数值频谱, 仅用于展示, 无实际用途
downSampledKspaceImage = rsos(downSampledKspace);  % RSOS 组合多个欠采样 k 空间, 重建得到 实值欠采样 k 空间频谱图像 (多对一) 128×120 double

figure;
imagesc(downSampledKspaceImage, [0, 10000]);
title('Downsampled k-space');
colormap('gray');
colorbar();

%% Check Downsampled Image
%%
%dsImage = ifft2(ifftshift(ifftshift(downSampledKspace, 2), 1));
dsBrainCoils = ifftshift(ifft2(fftshift(downSampledKspace)));  % 120×128×5 complex double
dsImage = rsos(dsBrainCoils);  % 120×128 double

figure;
imagesc(dsImage);
title('Downsampled Image');
colormap('gray');
colorbar();

%% Find True ACS
% Find the length of the true ACS domain, including lines due to downsampling.

%finder = any(downSampledKspace, [2, 3]);  % 原代码
%assert(isequal(finder, mean(mask, [2, 3])), 'The indices are wrong.')  % 原代码
finder = any(any(downSampledKspace, 2), 3);  % finder.size = 120×1 logical
assert(isequal(finder, mean(mean(mask, 2),3)), 'The indices are wrong.')  % mean.size = 120×1 double

acsFinder = zeros(phaseLength, 1);  % 120×1 double
for idx = 2:phaseLength-1
    if finder(idx-1) == 1 && finder(idx) == 1 && finder(idx+1) == 1  
        acsFinder(idx-1:idx+1) = 1;  % ACS 区域置 1
    end
end

% Getting the idices of the ACS lines.
acsLines = find(acsFinder);  % 13×1 double (54-66)
% 之所以为 13 而非 12 条, 是因为多出的 1 条来自奇数位置索引的已采样数据

% Checking whether the parameters fit.
matrixHeight = samplingRate .* (kernelHeight-1) + 1;  % 2 * (4-1) + 1 = 7  表示 kernel/权重覆盖范围的高度为 7
acsTrueLength = length(acsLines);  % 表示真实的 ACS 线长度为 13 

assert(acsTrueLength >= matrixHeight, 'ACS is too short for kernel params')
assert((acsLength <= acsTrueLength) && (acsTrueLength <= acsLength+2), 'Mask is incorrect')

%% Building the weight matrix
% Has equation of X * w = Y in mind. (Eq.11)
%     X = inMatrix, Y = outMatrix, w = weights.
% Each kernel adds one row to the X and Y matrices.

acsPhases = acsTrueLength - matrixHeight + 1;          % 13 - 7 + 1 = 7, kernel 纵向移动次数 / 起始行次序
numKernels = acsPhases .* freqLength;                  % 7 * 128 = 897, kernel 在 ACS 上移动的总次数 (训练对数)
kernelSize = kernelHeight .* kernelWidth .* numCoils;  % 4 * 3 * 5 = 60, kernel 包含的权重个数
outSize = numCoils .* (samplingRate - 1);              % 5 * (2 - 1) = 5, outMatrix 的列数 ???

% Eq.11 S represents the collected signal in each element coil at some position
inMatrix = zeros(numKernels, kernelSize);  % 896×60 complex double
% Eq.11 S(m) 
outMatrix = zeros(numKernels, outSize);    % 896×5 complex double

hkw = floor(kernelWidth/2);     % Half kernel width.  1
hkh = kernelHeight/2;           % Half kernel height. 2
kidx = 1;  % "Kernel index" for counting the number of kernels.

% acsline = 54,55,56,57,58,59,60, 表示 kernel 的起始行
for acsLine = acsLines(1:acsPhases)'  % 遍历第 1,2,3,4,5,6,7 位   (acsLines 整体范围从 54 到 66)
    % 第一个 phases = linsapce(54, 54+7-1, 4) = linspace(54, 60, 4) = 54  56  58  60
    phases = linspace(acsLine, acsLine+matrixHeight-1, kernelHeight);  % Phases of the kernel
    % 频率编码轴长度 freqLength = 128, freq 表示当前 kernel 中点的横坐标/频率值 (1-128)
    for freq = 1:freqLength
        % 第一个 freqs = linspace(1-1, 1+1, 3) = linspace(0, 2, 3) = 0  1  2
        freqs = linspace(freq-hkw, freq+hkw, kernelWidth);  % Frequencies of the kernel.
        % 第一个 freqs = mod((0  1  2) - 1, 128) + 1 = 128  1  2
        freqs = mod(freqs-1, freqLength) + 1;  % For circular indexing.
        
        % 第一个 selected = linspace(56+1, 58-1, 2-1) = 57
        selected = linspace(phases(hkh)+1, phases(hkh+1)-1, samplingRate-1);
        % 第一个 selected = mod(57-1, 120) + 1 = 57
        selected = mod(selected-1, phaseLength) + 1;  % Selected Y phases.
        
        % 第一个 tempX = downSampledKspace((54  56  58  60), (128  1  2), :) 
        % 共 4 * 3 * 5 = 60 个数据点
        tempX = downSampledKspace(phases, freqs, :);  % 4×3×5 complex double
        % 第一个 tempY = downSampledKspace(57, 1, :)
        % 共 1 * 1 * 5 = 5 个待拟合点
        tempY = downSampledKspace(selected, freq, :);  % 1×1×5 complex double
        
        % Filling in the matrices row by row.
        inMatrix(kidx, :) = reshape(tempX, 1, kernelSize);  % ACS 训练数据点排在第 kids 行上, size = 1×60
        outMatrix(kidx, :) = reshape(tempY, 1, outSize);    % ACS 待拟合数据点排在第 kids 行上, size = 1×5
        
        kidx = kidx + 1;  % 下一保存行, 共 7 * 128 = 896 行, 对应 896 个 ACS 待拟合数据点
    end
end

assert(isequal(size(inMatrix), [numKernels, kernelSize]), 'Incorrect X matrix');
assert(isequal(size(outMatrix), [numKernels, outSize]), 'Incorrect Y matrix');

% Has equation of X * w = Y in mind. (Eq.11)
%     X = inMatrix, Y = outMatrix, w = weights.
% Each kernel adds one row to the X and Y matrices.

% Eq.11 求解 n^(m) (weights) 
%     pinv(inMatrix)    60×896    complex double
%     outMatrix         896×5     complex double
weights = pinv(inMatrix) * outMatrix;  % Calculate the weight matrix.  60×5 complex double 
%     B = pinv(A) returns the Moore-Penrose pseudoinverse of A (伪逆)
                       
assert(isequal(size(weights), [kernelSize, outSize]), 'Incorrect weight size');

%% GRAPPA Reconstruction
% Performing a naive reconstruction according to first principles causes an 
% overlap problem.
% 
% The lines immediately before and after the ACS lines are not necessarily 
% spaced with the sampling rate as the spacing.
% 
% This causes alteration of the original data if GRAPPA reconstruction is 
% performed naively.
% 
% The solution is to perform reconstruction on a blank, and then overwrite 
% all vlaues with the original data.
% 
% This alleviates the problem of having to do special operations for the 
% values at the edges.
% 
% Also, the lines of k-space at the start or end of k-space may be neglected 
% (depending on the implementation).
% 
% This requires shifting the finder by the sampling rate to look at the phase 
% from one step above.
% 
% If the downsampling does not match due to incorrect k-space dimensions, 
% errors will be overwritten by the final correction process.

% 根据第一个原理进行朴素的重建会导致重叠问题. 
% ACS 线前后的线不一定要以采样率作为间隔来间隔。 
% 如果朴素地执行 GRAPPA 重建，将导致原始数据的更改。 
% 解决方案是对空白重建，然后用原始数据覆盖所有值。  (数据一致性)
% 这减轻了必须对边缘值进行特殊运算的问题。 
% 而且，可以忽略在 k 空间的开始或结尾处的 k 空间的线（取决于实现方式）。 
% 这需要将 finder 按采样率移动，以从上方一步查看相位。 
% 如果由于不正确的 k 空间尺寸导致下采样不匹配，则最终的校正过程将覆写错误。

% Find the indices to fill, including at the beginning of k-space.
temp1 = find(diff(finder) == -1) + 1;  % Gets nearly all the lines.
% 待填充行索引???      3、5 ... 51、53、67、69 ... 117、119 共 53 行
%   diff(finder) 在 54-65 行为 0, 其余奇数行为 1, 偶数行为 -1  119×1 double
%   find(diff(finder) == -1) + 1 找到值为 -1 的偶数行的索引, 然后将索引+1,  53×1 double

temp2 = find(diff(circshift(finder, samplingRate)) == -1) - samplingRate + 1;
% 待填充行索引???      1、3 ... 51、53、67、69 ... 115、117 共 53 行
%   circshift(finder, samplingRate) 偶数行及 56-68 行为 1, 其余行为 0, 120×1 logical
%   diff(circshift(finder, samplingRate)) 在 56-67 行为 0, 其余奇数行为 1, 偶数行为 -1, 119×1 double  
%   find(diff(circshift(finder, samplingRate)) == -1) - samplingRate + 1 找到值为 -1 的偶数行的索引, 
%      然后将索引 - 2 + 1 = -1, 53×1 double
                                                                   
% Second line exists to catch the first few empty lines, if any are present.
% 待重建从而填充的缺失行索引 (mask=0)
fillFinder = unique([temp1; temp2]);  % 1、3 ... 51、53、67、69 ... 117、119 共 54 行

% Shift from first fill line to beginning of kernel data.
upShift = (hkh-1) .* samplingRate + 1;  % (2-1) * 2 + 1 = 3
% Shift from first fill line to end of kernel data.
downShift = hkh .* samplingRate - 1;  % 2 * 2 - 1 = 3

% GRAPPA 重建 k 空间
grappaKspace = zeros(size(downSampledKspace));  % 120×128×5 complex double

% 遍历待重建填充的缺失行索引 (mask=0) 1×53 double
for phase = fillFinder'
    % 第一个 phases = linspace(1-3, 1+3, 4) = -2  0  2  4
    phases = linspace(phase-upShift, phase+downShift, kernelHeight);
    % 第一个 phases = mod((-2  0  2  4) - 1, 120) + 1 = 118  120  2  4
    phases = mod(phases-1, phaseLength) + 1;  % Circularly indexed phases.
    % 频率编码轴长度 freqLength = 128, freq 表示当前 kernel 中点的横坐标/频率值 (1-128)
    for freq = 1:freqLength    
        % 第一个 freqs = linspace(1-1, 1+1, 3) = linspace(0, 2, 3) = 0  1  2
        freqs = linspace(freq-hkw, freq+hkw, kernelWidth);  % Frequencies of the kernel.
        % 第一个 freqs = mod((0  1  2) - 1, 128) + 1 = 128  1  2
        freqs = mod(freqs-1, freqLength) + 1;  % Circularly indexed frequencies.   
        
        kernel = downSampledKspace(phases, freqs, :);  % 取出已知数据点, 4×3×5 complex double
        % One line of the input matrix.
        tempX = reshape(kernel, 1, kernelSize);  % 1×60 complex double
        
        % One line of the output matrix.  
        tempY = tempX * weights;   % 加权线性组合得到缺失点, 1×5 complex double
        tempY = reshape(tempY, (samplingRate-1), 1, numCoils);  % 1×1×5 complex double
        
        % Selected lines of the output matrix to be filled in.
        % 第一个 selected = linspace(0+1, 2-1, 2-1) = 1
        selected = linspace(phases(hkh)+1, phases(hkh+1)-1, samplingRate-1);
         % 第一个 selected = mod(1-1, 120) + 1 = 1
        selected = mod(selected-1, phaseLength) + 1;  % Selected Y phases.
        
        grappaKspace(selected, freq, :) = tempY;  % 缺失点估计值填充
    end
end

% Filling in all the original data.
% Doing it this way solves the edge overlap problem.
% 只估计未采集的数据, 而直接拷贝已采集的数据, 从而实现数据一致性 (data consistency)
grappaKspace(finder, :, :) = downSampledKspace(finder, :, :);  % 数据一致性处理

assert(isequal(size(downSampledKspace), size(grappaKspace)), 'Incorrect matrix size.');
% Somewhat redundant assertion now.
% assert(isequal(grappaKspace(finder, :, :), downSampledKspace(finder, :, :)));

%% Display recon image
%%
%recon = ifft2(ifftshift(ifftshift(grappaKspace, 2), 1));  源代码
recon = ifftshift(ifft2(fftshift(grappaKspace)));  % 120×128×5 complex double 修正
reconImage = rsos(recon);  % 120×128 double  RSOS 组合重建结果

imagesc(reconImage);
title('Reconstructed Image');
colormap('gray');
axis('off');
colorbar();

%% Display phase angle
%temp = recon(:, :, 1);  % 120×128 complex double
% angle(reconImage) 则由于 reconImage 是实数图像, 相位角已经为 0 了

% P = angle(Z) returns the phase angles, in radians, for each element of complex array Z.
% The angles lie between ±π.
% For complex Z, the magnitude R and phase angle theta are given by
%       R = abs(Z)
%       theta = angle(Z)
% and the statement
%       Z = R.*exp(i*theta)
% converts back to the original complex Z.

figure;
colormap('gray');
 
subplot(2, 3, 1);
imagesc(angle(recon(:, :, 1)));
title('coil-1 recon phase angle');
colorbar();
 
subplot(2, 3, 2);
imagesc(angle(recon(:, :, 2)));
title('coil-2 recon phase angle');
colorbar();
 
subplot(2, 3, 3);
imagesc(angle(recon(:, :, 3)));
title('coil-3 recon phase angle');
colorbar();
 
subplot(2, 3, 4);
imagesc(angle(recon(:, :, 4)));
title('coil-4 recon phase angle');
colorbar();

subplot(2, 3, 5);
imagesc(angle(recon(:, :, 5)));
title('coil-5 recon phase angle');
colorbar();
 
subplot(2, 3, 6);
imagesc(angle(reconImage));
title('RSOS recon MR phase angle');
colorbar();

%% Compute difference image
%%
deltaImage = reconImage - originalImage;  % difference image  重建图像与标签图像之差

figure;
imagesc(deltaImage);   
title('Difference Image');
colormap('gray');
axis('off');
colorbar();

%% Summary
%%
figure;

subplot(2, 2, 1);
imagesc(originalImage);
title('Original Image');
colormap('gray');
axis('off');
colorbar();

subplot(2, 2, 2);
imagesc(dsImage);
title('Down-Sampled Image');
colormap('gray');
axis('off');
colorbar();

subplot(2, 2, 3);
imagesc(reconImage);
title('Reconstructed Image');
colormap('gray');
axis('off');
colorbar();

subplot(2, 2, 4);
imagesc(deltaImage);
title('Difference Image');
colormap('gray');
axis('off');
colorbar();

toc

参考文献：https://zhuanlan.zhihu.com/p/240454014

参考链接：GitHub - veritas9872/Medical-Imaging-Tutorial: A few projects with MRI and X-ray CT imaging that I have done. SENSE, GRAPPA, and linear CT reconstruction are included.

何处闻韶

关注

13
点赞
踩
46

收藏

觉得还不错? 一键收藏
10
评论
【MRI】GRAPPA (GeneRalized Autocalibrating Partially Parallel Acquisitions) 算法仿真实验与原理剖析 (Matlab 实现)

参考文献https://zhuanlan.zhihu.com/p/240454014
复制链接

扫一扫