4.4.3 EPI伪影：模拟Nyquist ghost（matlab）

MRI 学习笔记，为了方便记录，序号和俎栋林老师、高家红老师编写的《核磁共振成像》一书的序号保持一致，方便翻书对照。
参考：
1、《Handbook of MRI Pulse Sequences》，MATT A. BERNSTEIN, KEVIN F. KING and XIAOHONG JOE ZHOU
2、《核磁共振成像——物理原理和方法》，俎栋林老师、高家红老师
3、叶恩茂. 医学磁共振EPI中Ghost伪影消除方法研究[D]. 第一军医大学, 2006.

EPI

平面回波成像序列： Echo Planar Imaging
EPI序列是在一次射频脉冲激励后，利用正反相位连续切换采集信号。严格意义上来讲，EPI不属于序列分类，算是一种信号采集技术，是一种信号读出模式。
磁共振弥散加权成像（Diffusion Weighted Imaging, DWI），是唯一能够活体检测水分子扩散情况的无创影像检查技术，而且也是临床应用最多的一个磁共振功能监测技术。

问题

Nyquist ghost有哪几种情况？
Nyquist ghost在图像域是什么样的？
K空间什么样的操作会造成Nyquist ghost？包括1/2、1/4Nyquist ghost？
背后原理是什么？

单激发和多激发

单激发 Single-shot	多激发 Multi-shot
O 分辨率受限	O 分辨率高
X 图像形变大	O 图像形变小
X 磁化率伪影重	O 磁化率伪影轻
O 对运动伪影不敏感	X 运动伪影重

单激发是一次射频脉冲激励以后，通过读出信号的方式采集完一副图像所需要的所有K空间相位编码线。最大的有点是快，可以冻结宏观的运动伪影。目前临床EPI序列应用最多的是做DWI弥散，采用的大部分是单激发。所以腹部弥散，可以完全自由呼吸做而运动伪影却不大，因为单激发对运动伪影不敏感。

EPI factor

EPI factor：确定每次激励时收集的Ky线的数量
EPI因子越大，SNR越低，同时还会影响磁化率伪影和图像变形，扫描仪产生的噪音
而EPI序列，如果回波链多了EPI factor也会导致EPI伪影和变形加大；另外，EPI序列没有采用多个180°重聚脉冲聚集信号，对磁化率敏感，容易产生磁化率伪影。

Nyquist ghost

由于磁场的不均匀性，以及EPI通过正反梯度切换快速的读出信号，产生的涡流和化学位移等因素，使得K空间的奇、偶回波之间容易积累相位差，在相位编码方向上产生Ghost伪影。导致图像的变形（distortion）。

奇偶回波产生的K空间偏移方向不同，从而出现相位编码方向得相位震荡，产生Ghost伪影

1、Single-shot 和 Multi-shot 的伪影区别

在 Single-shot 的情况下，图像产生1/2 Nyquist ghosts
在 Multi-shot 的情况下，图像可能会产生1/4 Nyquist ghosts，1/8 Nyquist ghosts，由采集方式决定

2、Three types of Nyquist ghosts

three types of Nyquist ghosts in single-shot EPI:
1）constant-phase ghost ：与空间位置无关的相位误差引起的伪影，比如零阶的涡流梯度场和偏离中心视野成像的频率失配
2）linear-phase ghost ：频率编码方向上的线性相位误差引起的伪影，比如频率编码方向上的线性涡流场，梯度群时延和梯度放大器的延迟
3）oblique-phase ghost：相位编码方向上的线性相位误差引起的伪影，通常发生在EPI斜采时，产生的原因可能有各方向梯度不一致的涡流特性、群时延以及交叉涡流项

1）constant-phase ghost 恒定相位鬼影

k-space data are shifted $\varphi$ （spatially constant）
$I^{\prime}(l, m)=I(l, m) \cos \varphi+i I\left(l, m-\frac{N_{y}}{2}\right) \sin \varphi$

2）linear-phase ghost 线性位移鬼影

k-space data are shifted $\delta k_{x} = ± u \Delta k_{x}$ （频率编码方向）
$I^{\prime}(l, m)=I(l, m) \cos \left(\frac{2 \pi u l}{N_{x}}\right)+i I\left(l, m-\frac{N_{y}}{2}\right) \sin \left(\frac{2 \pi u l}{N_{x}}\right)$
我们可以把K空间，沿着频率编码方向，奇数行向前平移一个点，偶数行向后平移一个点，来达到这种伪影效果。

3）oblique Nyquist ghost 斜奈奎斯特鬼影

k-space data are shifted $\delta k_{y} = ± v \Delta k_{y}$ （相位编码方向）
$I^{\prime}(l, m)=I(l, m) \cos \left(\frac{2 \pi v m}{N_{y}}\right)+i I\left(l, m-\frac{N_{y}}{2}\right) \sin \left(\frac{2 \pi v\left(m-N_{y} / 2\right)}{N_{y}}\right)$
我们可以把K空间，沿着相位编码方向，奇数行向前平移一个点，偶数行向后平移一个点，来达到这种伪影效果。

Because this ghost is typically observed in oblique EPI scans, it has been called the oblique Nyquist ghost.
Common sources for the oblique Nyquist ghosts include inconsistent eddy current characteristics and group delays among the physical gradient axes (also known as gradient anisotropy; Aldefeld and Bornert 1998) and cross-term eddy currents.

4) a mixture of the ghosts

The combination of a constant-phase ghost and a linear-phase ghost will shift the vertical nodal line away from the image center.
The combination of a linear-phase ghost and an oblique Nyquist ghost can produce a tilted nodal line.

Matlab模拟

Simulated Phantom

定义Read、Phase方向点数

clear;
Readout = 256;
Phaseout = 256;

生成头部模拟图，尺寸为128*128

[P0,~] = phantom('Modified Shepp-Logan',128);

重新构建图像尺寸，填零至256*256

P = zeros(256,256);
P((256-128)/2:(256-128)/2+127,(256-128)/2:(256-128)/2+127) = P0;

figure;% imagesc(abs(P))
subplot(1,3,1);imshow(abs(P),[]);title('Software head phantom')

生成K空间

kspace =fftshift(fft2(ifftshift(P)));
subplot(1,3,2);imshow(abs(kspace),[]);title('Simulated K-space')

从K空间恢复到原图

I0 = fftshift(ifft2(ifftshift(kspace)));
subplot(1,3,3);imshow(abs(I0),[]);title('Reconstructed image')

Image

1、constant-phase ghost

1）1/2 Nyquist Ghost

Constant = 0.2;
ImgConstantPhaseHalf = P;
for j = 1:Readout
    for i = 128:Phaseout-1
        ImgConstantPhaseHalf(i,j) = complex(P(i,j)*cos(Constant),P(i-Phaseout/2+1,j)*sin(Constant));
    end
end
for j = 1:Readout
    for i = 1:127
        ImgConstantPhaseHalf(i,j) = complex(P(i,j)*cos(Constant),P(i+Phaseout/2,j)*sin(Constant));
    end
end
% figure;imshow(abs(ImgConstantPhaseHalf),[]);
% title("1/2 Nyquist Ghost")

2）1/4 Nyquist Ghost

Constant = 0.2;
ImgConstantPhaseQuarter = P;
for j = 1:Readout
    for i = 128:Phaseout
        ImgConstantPhaseQuarter(i,j) = complex(P(i,j)*cos(Constant),P(i-Phaseout/4,j)*sin(Constant));
    end
end
for j = 1:Readout
    for i = 1:127
        ImgConstantPhaseQuarter(i,j) = complex(P(i,j)*cos(Constant),P(i+Phaseout/4,j)*sin(Constant));
    end
end
% figure;imshow(abs(ImgConstantPhaseQuarter),[]);
% title("1/4 Nyquist Ghost")
3）1/2+1/4 Nyquist Ghost
Constant = 0.2;
ImgConstantPhaseHalfQuarter = P;
for j = 1:Readout
    for i = 128:Phaseout-1
        ImgConstantPhaseHalfQuarter(i,j) = complex(P(i,j)*cos(Constant),P(i-Phaseout/4,j)*sin(Constant)+P(i-Phaseout/2+1,j)*sin(Constant));
    end
end
for j = 1:Readout
    for i = 1:127
        ImgConstantPhaseHalfQuarter(i,j) = complex(P(i,j)*cos(Constant),P(i+Phaseout/4,j)*sin(Constant)+P(i+Phaseout/2,j)*sin(Constant));
    end
end
% figure;imshow(abs(ImgConstantPhaseHalfQuarter),[]);
% title("1/2+1/4 Nyquist Ghost")

2、linear-phase ghost

LinearX = 1;
ImgLinearPhase = P;
for j = 1:Readout
    for i = 128:Phaseout-1
        ImgLinearPhase(i,j) = complex(P(i,j)*cos(2*pi*LinearX*j/Readout),P(i-Phaseout/2+1,j)*sin(2*pi*LinearX*j/Readout));
    end
end
for j = 1:Readout
    for i = 1:127
        ImgLinearPhase(i,j) = complex(P(i,j)*cos(2*pi*LinearX*j/Readout),P(i+Phaseout/2,j)*sin(2*pi*LinearX*j/Readout));
    end
end
% figure;imshow(abs(ImgLinearPhase),[]);
% title("Linear Phase")

3、oblique Nyquist ghost

LinearY = 1;
ImgObliquePhase = P;
for j = 1:Readout
    for i = 128:Phaseout-1
        ImgObliquePhase(i,j) = complex(P(i,j)*cos(2*pi*LinearY*i/Phaseout),P(i-Phaseout/2+1,j)*sin(2*pi*LinearY*(i-Phaseout/2)/Phaseout));
    end
end
for j = 1:Readout
    for i = 1:127
        ImgObliquePhase(i,j) = complex(P(i,j)*cos(2*pi*LinearY*i/Phaseout),P(i+Phaseout/2,j)*sin(2*pi*LinearY*(i-Phaseout/2)/Phaseout));
    end
end
% figure;imshow(abs(ImgObliquePhase),[]);
% title("Oblique Phase Ghost")

K-Space

1、constant-phase ghost

1）1/2 Nyquist Ghost

ShiftConstantPhase = zeros(Readout,Phaseout);
ShiftConstantPhase(1:2:256,:) = -0.2*1i;
ShiftConstantPhase(2:2:256,:) = 0.2*1i;
KspaceConstantHalf = kspace .* exp(ShiftConstantPhase);
% figure;
% imshow(abs(fftshift(ifft2(ifftshift(KspaceConstantHalf)))),[]);
% title('1/2 Nyquist Ghost')

2）1/4 Nyquist Ghost

ShiftConstantPhase = zeros(Readout,Phaseout);
ShiftConstantPhase(1:4:256,:) = -0.5*1i;
ShiftConstantPhase(2:4:256,:) = -0.5*1i;
ShiftConstantPhase(3:4:256,:) = 0.5*1i;
ShiftConstantPhase(4:4:256,:) = 0.5*1i;
KspaceConstantQuarter = kspace .* exp(ShiftConstantPhase);
% figure;
% imshow(abs(fftshift(ifft2(ifftshift(KspaceConstantQuarter)))),[]);
% title('1/4 Nyquist Ghost')

3）1/2+1/4 Nyquist Ghost

ShiftConstantPhase = zeros(Readout,Phaseout);
ShiftConstantPhase(1:4:256,:) = -0.5*1i;
ShiftConstantPhase(2:4:256,:) = 0.5*1i;
ShiftConstantPhase(3:4:256,:) = 0.5*1i;
ShiftConstantPhase(4:4:256,:) = 0.5*1i;
KspaceConstantHalfQuarter = kspace .* exp(ShiftConstantPhase);
% figure;
% imshow(abs(fftshift(ifft2(ifftshift(KspaceConstantHalfQuarter)))),[]);
% title('1/2+1/4 Nyquist Ghost')

4）Print Image

figure;
subplot(131);
imshow(abs(fftshift(ifft2(ifftshift(KspaceConstantHalf)))),[]);
title('1/2 Nyquist Ghost')
subplot(132);
imshow(abs(fftshift(ifft2(ifftshift(KspaceConstantQuarter)))),[]);
title('1/4 Nyquist Ghost')
subplot(133);
imshow(abs(fftshift(ifft2(ifftshift(KspaceConstantHalfQuarter)))),[]);
title('1/2+1/4 Nyquist Ghost')

2、linear-phase ghost

KspaceLinear = zeros(Phaseout,Readout);
for i = 2:Phaseout-1
    for j = 2:Readout-1
        if mod(i,2) == 0
            KspaceLinear(i,j) = kspace(i,j-1);
        else
            KspaceLinear(i,j) = kspace(i,j+1);
        end
    end        
end
ImgLinearPhase = fftshift(ifft2(ifftshift(KspaceLinear)));
figure;imshow(abs(ImgLinearPhase),[]);title('Linear Phase(X)')

3、oblique Nyquist ghost

KspaceOblique = zeros(Phaseout,Readout);
for i = 2:Phaseout-1
    for j = 1:Readout
        if mod(i,2) == 0
            KspaceOblique(i,j) = kspace(i-1,j);
        else
            KspaceOblique(i,j) = kspace(i+1,j);
        end
    end        
end
ImgObliquePhase = fftshift(ifft2(ifftshift(KspaceOblique)));
figure;imshow(abs(ImgObliquePhase),[]);title('Oblique Phase Ghost(Y)')

4、a mixture of the ghosts

KspaceMixture = zeros(Phaseout,Readout);
for i = 2:Phaseout-1
    for j = 1:Readout
        if mod(i,2) == 0
            KspaceMixture(i,j) = kspace(i-1,j)*exp(0.5*1i);
        else
            KspaceMixture(i,j) = kspace(i+1,j)*exp(-0.5*1i);
        end
    end        
end
ImgMixturePhase = fftshift(ifft2(ifftshift(KspaceMixture)));
% figure;imshow(abs(ImgMixturePhase),[]);title('Mixture Ghost')