11.2.1 重建圖像的復(fù)雜性

Physical spin density 純凈的位置圖像 p(x) 是一個嚴(yán)格的實數(shù)矩陣，但是重建的圖像 `p(x) 大多數(shù)時候是一個復(fù)數(shù) ，相位不是0。

最簡單的一個例子就是有一個全局相位偏移存在于空間頻域信號s(k)的每一個數(shù)據(jù)點上，導(dǎo)致原始的純凈的s(k) 被相位調(diào)制了，成為`s(k)：

全局相位的偏移，一般主要是由于實部通道和虛部通道被替換了，或者是不正確的反調(diào)制。這會直接導(dǎo)致在重建的圖像上也會有一個全局的位偏。 此處，s(k) 是 p(x) 的傅里葉變換得到的。

此時，僅僅依靠 `p(x) 的 實部 不能給出一個正確的重建圖像，除非我們可以理解那個相位偏移怎么產(chǎn)生的，并且修正這樣的偏移。

一般地，通常的做法是取`p(x)的幅度值，又叫做magnitude image，作為最后重建的圖像，這個圖像獨立于相位，可以很好地規(guī)避有許多來源的全局的相位偏移誤差的影響。

在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) 的取模結(jié)果。 這樣的local phase errors 可以被移除掉，像移除 global phase errors 一樣。這樣對重建的magnitude image 就沒有過多影響了。

我們也可以在k-space中引入一個線性位移，那么就變成了如下形式：

這樣在重建得到`p(x) 就會有一個空間位移x0,即 `p(x)---> `p(x+x0)

11.2.3 相位成像和相位混疊（包裹）

講完圖像重建和位移定理，下面我們來講相位成像（phase imaging）。 盡管在MRI 成像中，經(jīng)常使用magnitude image 來表示2D MRI 重建圖像，但是我們丟棄的相位信息有時候是需要使用的。因此，畫一個2D的phase image是有必要的， 它的pixel intensity 是正比于計算得到的complex `p(x)的相位值的。為了獲得這樣的一個相位圖，我們經(jīng)常使用這樣一個計算式：

反正切函數(shù)可以用來求2D pixel intensity值，但是有一個問題，反正切函數(shù)是一個在[-pi,pi)上的周期函數(shù)，所有的相位值都會被映射到這個區(qū)間，盡管真實的相位值，可以取任何的有限實數(shù)值。這就會導(dǎo)致，相位值相差多個2pi的相位，在[-pi,pi)這個區(qū)間里面被壓縮或者相同表示，導(dǎo)致部分相位失序，從而造成相位混疊（phase aliasing）。常見的相位混疊是帶狀或者是斑馬條紋那樣的偽跡。

再次考察 s(k)-----> s(k-k0)，至于這樣的位移發(fā)生的原因，后面會提到，這里就按下不表。信號在空間頻域有一個k0的偏移，導(dǎo)致：

如此，多了一個相位，該相位隨著k0和x變化而變化。導(dǎo)致一些相位值超出[-pi,pi)的范圍，最后被強(qiáng)制映射回這個區(qū)間，造成混疊（相位包裹）。k0越大，生成的混疊帶狀就越多。a)圖是原始的幅度譜，b)圖是原始的相位譜， c) 圖是原始s(k)偏移一個單位得到的相位混疊圖，? d)圖是原始s(k)偏移5個單位得到的相位混疊圖。

11.2.4 對偶性

這里就是2個簡單的傅里葉變換對：

偏移的delta 函數(shù)，有一個周期變化，當(dāng)k--->k+1/x0? H(k) 不變，這個與數(shù)據(jù)采樣會有關(guān)系，后面會提到。

最后的最后，感謝大家看完枯燥的MRI 傅里葉圖像變化，不過應(yīng)該相比于信號處理那本書上要更加形象，畢竟有圖有真相，哈哈哈，關(guān)鍵是看到這些變換的具體應(yīng)用，心里更加踏實，而不是停留在simulation的基礎(chǔ)之上。

這里大家可以玩一個小實驗，把兩個圖的magnitude image 和phase image 調(diào)換，看看 最后生成什么結(jié)果。

大家可以玩一下，蠻好玩的，哈哈哈，代碼無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)')