影象和多維 FT
在醫學成像,光譜學,影象處理,密碼學和其他科學和工程領域中,通常情況是希望計算影象的多維傅立葉變換。這在 Matlab 中非常簡單:(多維)影象只是 n 維矩陣,畢竟傅立葉變換是線性運算元:一個只是迭代傅立葉變換其他維度。Matlab 提供了 fft2
和 ifft2
來實現 2-d 或 fftn
的 n 維。
一個潛在的缺陷是影象的傅立葉變換通常顯示為中心有序,即影象中間的 k 空間的原點。Matlab 提供 fftshift
命令以適當地交換傅立葉變換的 DC 分量的位置。這種排序符號使得執行常見的影象處理技術變得非常容易,其中之一如下所示。
零填充
將小影象內插到更大尺寸的一種快速且骯髒的方法是對其進行傅立葉變換,用零填充傅立葉變換,然後進行逆變換。這有效地在具有 sinc 形基函式的每個畫素之間插值,並且通常用於向上擴充套件低解析度醫學成像資料。讓我們從載入內建影象示例開始
%Load example image
I=imread('coins.png'); %Load example data -- coins.png is builtin to Matlab
I=double(I); %Convert to double precision -- imread returns integers
imageSize = size(I); % I is a 246 x 300 2D image
%Display it
imagesc(I); colormap gray; axis equal;
%imagesc displays images scaled to maximum intensity
我們現在可以獲得 I 的傅立葉變換。為了說明 fftshift
的作用,讓我們比較兩種方法:
% Fourier transform
%Obtain the centric- and non-centric ordered Fourier transform of I
k=fftshift(fft2(fftshift(I)));
kwrong=fft2(I);
%Just for the sake of comparison, show the magnitude of both transforms:
figure; subplot(2,1,1);
imagesc(abs(k),[0 1e4]); colormap gray; axis equal;
subplot(2,1,2);
imagesc(abs(kwrong),[0 1e4]); colormap gray; axis equal;
%(The second argument to imagesc sets the colour axis to make the difference clear).
我們現在已經獲得了示例影象的 2D FT。為了對它進行零填充,我們想要獲取每個 k 空間,用零填充邊緣,然後進行後向變換:
%Zero fill
kzf = zeros(imageSize .* 2); %Generate a 492x600 empty array to put the result in
kzf(end/4:3*end/4-1,end/4:3*end/4-1) = k; %Put k in the middle
kzfwrong = zeros(imageSize .* 2); %Generate a 492x600 empty array to put the result in
kzfwrong(end/4:3*end/4-1,end/4:3*end/4-1) = kwrong; %Put k in the middle
%Show the differences again
%Just for the sake of comparison, show the magnitude of both transforms:
figure; subplot(2,1,1);
imagesc(abs(kzf),[0 1e4]); colormap gray; axis equal;
subplot(2,1,2);
imagesc(abs(kzfwrong),[0 1e4]); colormap gray; axis equal;
%(The second argument to imagesc sets the colour axis to make the difference clear).
此時,結果相當不起眼:
然後我們進行反向變換後,我們可以看到(正確!)零填充資料為插值提供了一種合理的方法:
% Take the back transform and view
Izf = fftshift(ifft2(ifftshift(kzf)));
Izfwrong = ifft2(kzfwrong);
figure; subplot(1,3,1);
imagesc(abs(Izf)); colormap gray; axis equal;
title('Zero-filled image');
subplot(1,3,2);
imagesc(abs(Izfwrong)); colormap gray; axis equal;
title('Incorrectly zero-filled image');
subplot(1,3,3);
imagesc(I); colormap gray; axis equal;
title('Original image');
set(gcf,'color','w');
請注意,零填充影象大小是原始影象大小的兩倍。在每個維度中,可以將填充零填充超過兩倍,但顯然這樣做不會隨意增加影象的大小。
提示,技巧,3D 及以上
上述示例適用於 3D 影象(例如,通常由醫學成像技術或共聚焦顯微鏡產生),但是例如需要 fft2
替換為 fftn(I, 3)
。由於多次編寫 fftshift(fft(fftshift(...
有點繁瑣,因此在本地定義 fft2c
等函式以在本地提供更簡單的語法是很常見的 - 例如:
function y = fft2c(x)
y = fftshift(fft2(fftshift(x)));
請注意,FFT 很快,但是大型的多維傅立葉變換仍然需要在現代計算機上進行。它還具有內在的複雜性:上面顯示了 k 空間的大小,但相位絕對至關重要; 影象域中的平移等效於傅立葉域中的相位斜坡。人們可能希望在傅立葉域中進行一些更復雜的操作,例如過濾高或低空間頻率(通過將其與濾波器相乘),或者遮蔽掉與噪聲相對應的離散點。相應地,大量社群生成的程式碼用於處理 Matlab 的主社群儲存庫站點 File Exchange 上可用的常見傅立葉操作。