11.2.1 重建图像的复杂性
Physical spin density 纯净的位置图像 p(x) 是一个严格的实数矩阵，但是重建的图像 `p(x) 大多数时候是一个复数 ，相位不是0。
最简单的一个例子就是有一个全局相位偏移存在于空间频域信号s(k)的每一个数据点上，导致原始的纯净的s(k) 被相位调制了，成为`s(k)：
全局相位的偏移，一般主要是由于实部通道和虚部通道被替换了，或者是不正确的反调制。这会直接导致在重建的图像上也会有一个全局的位偏。 此处，s(k) 是 p(x) 的傅里叶变换得到的。
此时，仅仅依靠 `p(x) 的 实部 不能给出一个正确的重建图像，除非我们可以理解那个相位偏移怎么产生的，并且修正这样的偏移。
一般地，通常的做法是取`p(x)的幅度值，又叫做magnitude image，作为最后重建的图像，这个图像独立于相位，可以很好地规避有许多来源的全局的相位偏移误差的影响。
在MRI 图像重建的最后，一般都是取magnitude image， 这在惯常的操作，但是，对于local phase error 我们需要采取其他的方法 来修正。
11.2.2 偏移定理
考虑一个信号 s(k) 有一个位置偏移k0， 使得 s(k) ----> s(k-k0) 。这里信号s(k)的中心就从移动到了k=k0位置。咱们看一下，这个位置偏移对于 `p(x) 图像重建会有什么影响。对sm(k-k0) 施加反傅里叶变换，得到估计的`p(x) 图像。
和原来`p[expected](x) 相比起来，多了一个相位，这样的话，如果对 `p(x) 取模运算，得到的magnitude image 就还是原来的`p[expected](x) 的取模结果。 这样的local phase errors 可以被移除掉，像移除 global phase errors 一样。这样对重建的magnitude image 就没有过多影响了。
我们也可以在k-space中引入一个线性位移，那么就变成了如下形式：
这样在重建得到`p(x) 就会有一个空间位移x0,即 `p(x)---> `p(x+x0)
11.2.3 相位成像和相位混叠（包裹）
讲完图像重建和位移定理，下面我们来讲相位成像（phase imaging）。 尽管在MRI 成像中，经常使用magnitude image 来表示2D MRI 重建图像，但是我们丢弃的相位信息有时候是需要使用的。因此，画一个2D的phase image是有必要的， 它的pixel intensity 是正比于计算得到的complex `p(x)的相位值的。为了获得这样的一个相位图，我们经常使用这样一个计算式：
反正切函数可以用来求2D pixel intensity值，但是有一个问题，反正切函数是一个在[-pi,pi)上的周期函数，所有的相位值都会被映射到这个区间，尽管真实的相位值，可以取任何的有限实数值。这就会导致，相位值相差多个2pi的相位，在[-pi,pi)这个区间里面被压缩或者相同表示，导致部分相位失序，从而造成相位混叠（phase aliasing）。常见的相位混叠是带状或者是斑马条纹那样的伪迹。
再次考察 s(k)-----> s(k-k0)，至于这样的位移发生的原因，后面会提到，这里就按下不表。信号在空间频域有一个k0的偏移，导致：
如此，多了一个相位，该相位随着k0和x变化而变化。导致一些相位值超出[-pi,pi)的范围，最后被强制映射回这个区间，造成混叠（相位包裹）。k0越大，生成的混叠带状就越多。a)图是原始的幅度谱，b)图是原始的相位谱， c) 图是原始s(k)偏移一个单位得到的相位混叠图， d)图是原始s(k)偏移5个单位得到的相位混叠图。
11.2.4 对偶性
这里就是2个简单的傅里叶变换对：
偏移的delta 函数，有一个周期变化，当k--->k+1/x0 H(k) 不变，这个与数据采样会有关系，后面会提到。
最后的最后，感谢大家看完枯燥的MRI 傅里叶图像变化，不过应该相比于信号处理那本书上要更加形象，毕竟有图有真相，哈哈哈，关键是看到这些变换的具体应用，心里更加踏实，而不是停留在simulation的基础之上。
这里大家可以玩一个小实验，把两个图的magnitude image 和phase image 调换，看看 最后生成什么结果。
大家可以玩一下，蛮好玩的，哈哈哈，代码无bug，只需要替换图像即可。
clc;clear;
%% read the image
A=imread('C:\Users\LENOVO\Desktop\menalisa.jpeg');
B=imread('C:\Users\LENOVO\Desktop\wangdeming.jpeg');
% convert to the gray image
A_gray=double(rgb2gray(A));
B_gray=double(rgb2gray(B));
% figure()
% subplot(121)
% imshow(A_gray,[])
% title('Monalisa')
% subplot(122)
% imshow(B_gray,[])
% title('Wangdeming')
%% resize the image to 500*500
A_resize=imresize(A_gray,[400 400]);
B_resize=imresize(B_gray,[400 400]);
figure(1)
subplot(121)
imshow(A_resize,[])
title('A')
subplot(122)
imshow(B_resize,[])
title('B')
% Fourier trnasform on A and B
fA=fft2(A_resize);
fB=fft2(B_resize);
% get the phase and magnitude mapping
fA_mag=abs(fftshift(fA)); % center the low frequency area
fA_phase=angle(fA);
fB_mag=abs(fftshift(fB)); % center the low frequency area
fB_phase=angle(fB);
% scale the magnitude image
S_A=log(1+abs(fA_mag));
S_B=log(1+abs(fB_mag));
% only for visualization , apply normalization on the data
figure(2)
subplot(221)
imshow(S_A,[])
title('A magnitude image')
subplot(222)
imshow(fA_phase,[])
title('A phase image')
subplot(223)
imshow(S_B,[])
title('B magnitude image')
subplot(224)
imshow(fB_phase,[])
title('B phase image')
%% exchange the mag and phase
fA_mag_ifftshift=ifftshift(fA_mag);
fB_mag_ifftshift=ifftshift(fB_mag);
fA_cross=fA_mag_ifftshift.*cos(fB_phase)+fA_mag_ifftshift.*sin(fB_phase).*1i;
fB_cross=fB_mag_ifftshift.*cos(fA_phase)+fB_mag_ifftshift.*sin(fA_phase).*1i;
fA=uint8(abs(ifft2(fA_cross)));
fB=uint8(abs(ifft2(fB_cross)));
figure(3)
subplot(121)
imshow(fA,[])
title('reconstruction from mag(A) & phase(B)')
subplot(122)
imshow(fB,[])
title('reconstruction from mag(B) & phase(A)')
网友评论