轉自：??????CT圖像重構方法詳解——傅里葉逆變換法、直接反投影法、濾波反投影法_Absolute Zero-CSDN博客_反投影法

緒

? ? ? 在做CT圖像處理的時候遇到很多問題，對于濾波反變換有許多細節存在疑問，經過多天查找資料和利用MATLAB程序一步步實現后終于豁然開朗，于是想要總結成文，作為筆記方便今后查看。文中若有錯誤歡迎指出！

目 錄

（1）傅里葉逆變換法

（2）直接反投影法

（3）濾波反投影法

（4）MATLAB程序

---

? ? ? CT圖像的形成和重構，在數學上的描述分別為拉冬變換和拉冬反變換(RANDON:氡）。簡單點說，拉動變換是用于描述通過X射線掃描物體形成CT圖像這個過程的一種數學表示，而拉東反變換描述的則是對CT投影數據進行重構，還原成物體圖像的一種數學方法。本文章主要講的是CT圖像重構的方法，也就是拉東反變換的具體實現方法，當然可能也有小伙伴想了解一下什么是拉東變換，這里貼上傳送門：拉冬變換。(Radon（拉冬）變換可以用來進行直線識別的算法。它是一種原始灰度圖像到（p,o）二維矩陣的映射，映射關系是灰度圖像的像素值對參數（p,o）的直線的線積分，或者叫投影，可理解為垂直于這個直線方向的像素值累計求和)

其中拉冬反變換在數學上的表達式如下：

由公式可看出拉冬逆變換在數學上實際上是一個二維傅里葉逆變換，但由于二維傅里葉逆變換計算量大，耗時長，因此計算機對于拉冬逆變換的具體實現方法主要為反投影法，反投影法又分為直接反投影法和濾波反投影法。由此可見CT重構方法主要為以下三種：

（1）傅里葉逆變換法

（2）直接反投影法

（3）濾波反投影法

下面將具體對這三種圖像重構方法進行解釋。

---

（1）傅里葉逆變換法

該方法基于一個重要的定理：中心切片定理。

該定理簡單地理解就是：通過角度為θ掃描得到的投影，該投影的一維傅里葉變換，與對整個圖像二維傅里葉變換后，二維頻域中對應θ角度的一個切片信號是相同的，下面兩個圖理解起來更直觀。

圖1

圖2

根據該理論，傅里葉逆變換法可以簡單分成以下步驟：

① 假設每旋轉1°就掃描一次，當對物體掃描了180°之后，我們就能得到180個投影信號（就是180根投影線）→在臨床上，若使用平行掃描CT，我們拿到手的數據就是這個（在數學上，就是對圖像進行拉冬變換）

② 對180個投影信號進行一維傅里葉變換

③ 對②得到的180個一維頻域信號，根據相對應的掃描角度，在空間中旋轉排列，拼成一個二維頻域空間（如圖3）

④ 由于數據是離散的，直接按照角度進行排列難以鋪滿整個二維空間，因此還需要對空缺的地方進行插值（一般三次樣條插值效果最好），但插值會帶來一定誤差。另外，由于中心的信號密集，周圍的信號稀疏，顯然會損失一部分高頻數據，造成高頻信號失真，這就是采用傅里葉逆變換法重構圖像時會使得圖像邊緣模糊的原因。

⑤ 對④中拼接而成的二維圖像進行二維傅里葉逆變換，就可重構原圖

圖3? 傅里葉逆變換法

該方法的缺點有：

a/高頻信號有所失真

b/在插值時還涉及到極坐標和直角坐標的變換，計算量大

c/需要用到二維傅里葉逆變換，總體耗費時間長

由于計算機處理二維傅里葉逆變換的計算量太大，因此很少直接使用該方法實現拉東逆變換。

（2）直接反投影法

反投影法的原理是將所測得的投影值，按照其原投影路徑，平均地分配到經過的每一個點上，把各個方向的投影值都這樣反投影后，在把每個角度的反投影圖像進行累加，從而推斷出原圖。為方便理解，下面分別用兩張圖，通過MATLAB編程來還原一下該過程（程序在最后）：

原圖如下：（圖1簡單，圖2復雜）

①首先看一下圖像經過拉東變換（即CT掃描）后的數據（圖中每一列代表對應每一個角度的投影信號，在這里只是拼起來成一張圖了）

②對每個角度的信號進行反投影（選了幾個代表性角度）：

② 把以上角度發反投影圖像疊加起來：

當越來越多的反投影圖像加起來之后，可以看到：

到此，直接反投影法結束。值得注意的是，由于反投影圖是離散疊加的，顯然在中心處信號集中，邊緣處信號稀疏，因此在最后需要在空缺的地方進行插值，才能得到最終的原圖像。上面演示的過程中沒有進行插值，僅僅是把反投影圖疊加，因此能夠看出每一張圖都不是很光滑。

從以上過程可以很直觀地看出直接反投影法出現偽跡的原因：

①在“回抹“過程中（就是把投影信號平均到每個二維空間點的過程中），會把原圖像本來是0的像素點也”抹“上一個平均值，最終使得重構的圖像中存在誤差；

②插值過程中會帶來誤差，且周圍信號稀疏，高頻信號有所失真，導致圖像邊緣模糊；

③在反投影圖不斷疊加過程中，能夠看到這種疊加方式會帶來明顯的“星狀偽跡“，這是造成圖像邊緣模糊的最重要因素。在投影數據少的時候更明顯（現實中投影數據的多少取決于機器）。

關于第③點，在數學上有研究表明，反投影重構后的圖像fb(x,y)，和真實原圖f(x,y)之間存在以下卷積關系：

如果要重構成真實的圖像，就需要把這個由于反投影法本身帶來的1/r效應去掉，由于1/r會使得圖片變得模糊，因此1/r也稱之為模糊因子。在實際操作中，把1/r影響去除的方法這就是下面這種最常用的（計算機、商業普遍使用）CT圖像重構方法：濾波反投影法。

（3）濾波反投影法

從上面的式（1）可以看到，1/r跟原圖像是卷積關系，通過傅里葉變換轉變到頻域上后會變成乘積關系，這說明直接在頻率上處理會更加簡便。

在頻域中，式（3）中的r稱之為權重因子，其實就相當于一個濾波器，想要去除模糊因子的影響，重構成原始圖像，只需要三步：

①把重構圖像fb(x,y)進行二維傅里葉變換得到Fb(x,y)

②在頻域中把Fb(x,y)與濾波器r相乘

③把r×Fb(x,y)進行二維傅里葉逆變換即能夠得到原圖f(x,y)

具體過程見下圖，圖中的ρ就是公式中的r。

? ? ??

圖4? 直接反投影法

這種先反投影，后濾波的方法需要用到兩次二維傅里葉變換，計算量太大，因此也不是一個特別理想的矯正方法。

而濾波反投影法，顧名思義就是先濾波，后反投影的重構方法。具體步驟如下：

①對180個投影信號先分別進行一維傅里葉變換

②在頻域中對所有投影信號進行濾波，也就是乘上權重因子r

③把所有濾波后的信號再進行一維傅里葉逆變換，還原到時域

④對每一個已經濾波的投影信號進行反投影，最后疊加（這一步跟前面的直接反投影法步驟完全相同，只是投影信號已經過了濾波）

具體流程見下圖：

圖5? 濾波反投影法

這個方法的好處是，把兩次二維傅里葉變換變成了兩次一維傅里葉變換，計算速度大大提升。于是濾波反投影法的核心問題就變成了如何選擇一個合適的濾波器r。數學中的濾波器r是一個理想化的濾波函數，現實中不存在，因此只能設計一個近似的濾波函數來代替r的作用。

下面常用的濾波函數有：R-L濾波函數和S-L濾波函數。

①R-L濾波函數（Ram-Lak濾波器 Ramp Filter）

從頻域圖中可以看到該濾波器的作用是把低頻信號減少，從而突出高頻信號，因此經過該函數濾波后，重建圖像的輪廓會更清晰，并且函數簡單，實現起來更方便。但由于該函數在頻域上有一個加窗處理（就是把高頻部分直接截斷了），因此在時域中會出現有許多振蕩的小尾巴（截窗振蕩效應），這將會讓重構出來的圖像存在Gibb’s現象（重構圖像存在振蕩，不連續）。

②S-L濾波函數（Shepp-Logan濾波器）

S-L濾波器在頻域上不是直接加窗截斷，而是通過一些比較平滑的窗函數對函數進行約束。因此經過S-L濾波后重構的圖像振蕩更少，重構的圖像質量比R-L濾波器更好一些，但是因為S-L在高頻段并沒有直接截斷，偏離了理想濾波器的效果，因此在高頻段（輪廓）上的重構效果沒有R-L濾波器好。

可以在下圖更直觀地看到兩個濾波器之間的區別：

下面通過MATLAB編程來還原一下濾波反投影法的過程（程序在最后，只用復雜那個圖進行展示）：

①對每一列投影信號分別在頻域進行R-L濾波(高頻信號明顯突出了)

②對每個濾波后的投影信號反投影（這里開始與直接反投影法一樣）

③ 把所有反投影圖像疊加起來?

為了驗證手動重構的效果，下面將使用matlab中的iradon函數直接對CT投影數據進行重構，把手動重構的圖像和matlab自帶函數重構的圖像進行對比。

（iradon是matlab中的拉東逆變換函數，但是其實際采用的重構方法并不是對數據進行二維傅里葉逆變換，而是通過濾波反投影法實現的，且默認插值方法為“線性插值linear“，默認濾波方法為”R-L濾波函數“，這在matlab的help中可以直接看到。）

仔細觀察兩個重構結果，會發現iradon的重構結果更光滑一些，那是因為手動方法只是簡單粗暴地把全部離散的反變換圖片疊加起來，而iradon則會在最后對圖像進行線性插值，使得圖像更連續更光滑，重構效果更好。

最后把iradon重構的圖像放大來看看圖像細節：

可以看到重構的圖像其實并不連續，甚至能夠看出一些振蕩的波紋，這就是R-L濾波函數截窗振蕩效應所帶來的Gibb’s現象。

（4）MATLAB程序

%% CT反投影法 %%==========================Part 1========================== %% 1、直接反投影法 % 繪制一張簡單的圖，進行直接反投影法 clc,clear,close all % ================================================ % 簡單圖 % I=zeros(512,512); % for i=1:512 % for j=1:512 % if ((i-256)*(i-256)+(j-256)*(j-256))<1024 % I(i,j)=1; % end % end % end % figure,imshow(I,[]); % ===============================================%============================================ % 復雜圖 I = phantom(512); figure,imshow(I,[]); %============================================% 繪制代表角度的反投影圖像 R=radon(I,0:179); %對物體CT掃描180°的數據，每1°掃描一次 figure,imshow(R,[]);title('CT圖像'); % 繪制1、45、90、135、180的投影信號和反投影圖像 theta=[1 45 90 135]; for i=1:length(theta)r=R(:,i);II=iradon([r r],[theta(i) theta(i)],'linear','none')/2;%回抹（反投影）之前不濾波II1{i}=II;figure,subplot(1,2,1),imshow(r,[]);title([num2str(theta(i)) '°投影信號']);subplot(1,2,2),imshow(II,[]);title(['往' num2str(theta(i)) '°方向反投影']); end II2=II1{1}+II1{2}+II1{3}+II1{4}; figure,imshow(II2,[]);% 對反投影圖像進行疊加 r=R(:,1); II_r=iradon([r r],[1 1])/2; k=[30 15 10 1]; for j=1:4A=II_r;for i=2:k(j):180r=R(:,i);II=iradon([r r],[i i],'linear','none')/2;%回抹（反投影）之前不濾波A=A+II; end % figure,imshow(A,[min(min(I)) max(max(I))]);figure,imshow((A),[]); end%% =========================Part 2============================= %% 2、濾波反投影法 clc,clear,close all % 復雜圖 I = phantom(512); figure,imshow(I,[]);% 對投影做傅里葉變換 R=radon(I,0:179); width=length(R); % 設計R-L濾波器 filter=2*[0:round(width/2-1), width/2:-1:1]'/width; % 展示濾波前的CT投影信號 figure,imshow(R,[]),title('濾波前的CT投影信號');% 每一列做傅里葉變換 r_fft=fft(R,729); % 每一列做傅里葉變換后濾波 r_fft_filter=r_fft.*filter; % 濾波后反變換(real取實部) r_fft_filter_v=real(ifft(r_fft_filter)); % 展示濾波后的CT投影信號 % figure,imshow(r_fft_filter_v,[]),title('濾波后的CT投影信號'); figure,imshow(r_fft_filter_v),title('濾波后的CT投影信號');% 繪制1、45、90、135、180的投影信號和反投影圖像 theta=[1 45 90 135]; for i=1:length(theta)r=r_fft_filter_v(:,i);II=iradon([r r],[theta(i) theta(i)],'linear','none')/2;%由于iradon默認有濾波，因此關掉默認的濾波選項（none），直接用上面已經手動濾波的數據進行反投影II1{i}=II;figure,subplot(1,2,1),imshow(r,[]);title([num2str(theta(i)) '°投影信號']);subplot(1,2,2),imshow(II,[]);title(['往' num2str(theta(i)) '°方向反投影']); end II2=II1{1}+II1{2}+II1{3}+II1{4}; figure,imshow(II2,[]);% 對反投影圖像進行疊加 r=r_fft_filter_v(:,1); II_r=iradon([r r],[1 1])/2; k=[30 15 10 1]; for j=1:4A=II_r;for i=2:k(j):180r=r_fft_filter_v(:,i);II=iradon([r r],[i i],'linear','none')/2;%由于iradon默認有濾波，因此關掉默認的濾波選項（none），直接用上面已經手動濾波的數據進行反投影A=A+II; end % figure,imshow(A,[min(min(I)) max(max(I))]);figure,imshow((A),[]); end% 上面是反投影疊加的分步計算，下面直接輸入全部投影數據到iradon函數中進行驗證 B=iradon(r_fft_filter_v,0:179,'linear','none'); figure, imshow(B,[]);title('驗證：濾波后圖像');

---

以上就是對CT圖像重構的整理結果，更多是自己在學習過程中的理解，若有表述不當的地方歡迎指出。

-------------------------------------------------------?2021.05.27? 筆記補充 -----------------------------------------------------------

重構圖像（514×514）與原圖像（512×512）尺寸不一致的原因分析：

首先需要感謝樓下?study_art?同學提出了這個問題，以及?xjch1900?同學的提醒，一開始確實沒考慮過這個問題，后來查看源代碼找到了原因，在此更新一下筆記。

①radon函數的matlab文件解釋：

% Grandfathered syntax % R = RADON(I,THETA,N) returns a Radon transform with the % projection computed at N points. R has N rows. If you do not % specify N, the number of points the projection is computed at % is: % % 2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3 % % This number is sufficient to compute the projection at unit % intervals, even along the diagonal.

我們一開始先使用了radon函數來模擬了CT的掃描過程，得到的結果是每個角度的投影數據。假設我們掃描的是一張正方形的圖像，大小為512×512，如果我們要把掃描的數據都存儲下來，就必須設置一個足夠長的數組來記錄這些數據，顯然，當掃描方向為45°角時得到的掃描數據是最長的（如下圖所示），需要的數組長度至少為?512×√2≈725?。

回到matlab的radon函數，從上面的文件解釋中可看到radon的語法中的N便是用來定義這個最長長度的。但當我們沒有對這個N有專門的限制時，該函數就會自動計算出這個“最長長度”，為了涵蓋圖像形狀的不同情況（如有的圖像為長方形），它設置了一個公用的計算公式：

N=2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3

【其中ceil函數表示“向上取整”；floor表示“向下取整”；norm表示“范數計算”，用在這里其實就是計算正方形斜邊的長度】

如果套入這個公式，計算出來的N為：

size(I)=[512 512]; N=2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3N =729

這比自己的理論值725大了一些，但實際上這是為了涵蓋所有圖像情況造成的，這并不會對最后的結果帶來太大影響。

在正文最后的matlab程序中，我并沒有專門定義這個N，因此radon函數自動把N通過公式計算出來了，大小為729，同時還定義了180個掃描角度，所以得到的所有投影數據的矩陣大小為729×180.

②iradon函數的matlab文件解釋：

% I = IRADON(R,THETA,INTERPOLATION,FILTER,FREQUENCY_SCALING,OUTPUT_SIZE) % OUTPUT_SIZE is a scalar that specifies the number of rows and columns in % the reconstructed image. If OUTPUT_SIZE is not specified, the size is % determined from the length of the projections: % % OUTPUT_SIZE = 2*floor(size(R,1)/(2*sqrt(2))) % % If you specify OUTPUT_SIZE, IRADON reconstructs a smaller or larger % portion of the image, but does not change the scaling of the data.

在使用iradon函數進行圖像重構時，同樣有一個output_size的參數可以設置（表示輸出圖像的大小），我們可以直接設置為原始圖像的大小（即512），這樣得到的圖像就跟原圖像大小一模一樣了。但在正文最后部分，我同樣沒有專門設置圖像大小，所以該函數就會自動套入一條涵蓋了所有圖像形狀的計算公式：

OUTPUT_SIZE = 2*floor(size(R,1)/(2*sqrt(2)))

由于我們在第一步用radon函數模擬CT掃描時沒有專門定義N，得到的投影數據的矩陣R大小為729×180，因此直接套入這個output_size公式后就直接算出：

size(R,1)=729; OUTPUT_SIZE = 2*floor(size(R,1)/(2*sqrt(2)))OUTPUT_SIZE =514

所以最后得到的重構圖像大小就變成了514×514.

③解決方法

• ?使用radon函數時，自己定義N。如在正文的例子中，直接在radon函數中輸入N的理論計算值725，這樣在重構iradon時就不需要設置output_size了，輸出結果直接為512×512
• 使用radon函數時，不定義N，讓函數自己計算圖像的N，然后在重構iradon時設置output_size的大小為原始圖像大小512，輸出結果同樣為512×512

?備注：若radon和iradon函數都不專門設置時，會出現重構圖像與原圖像尺寸大小不一致的情況，但這并不會改變重構圖像本身的縮放比例，因此若沒有專門的尺寸大小必須前后一致的嚴格要求，可以不用專門設置。

總結

以上是生活随笔為你收集整理的【转】CT图像重构方法详解——傅里叶逆变换法、直接反投影法、滤波反投影法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。