《MATLAB應用圖像處理》課件第3章

3.1傅立葉變換3.2離散余弦變換(DCT)3.3射線投影積分變換習題

第3章MATLAB圖像變換

3.1傅?立?葉?變?換

3.1.1連續(xù)傅立葉變換信號在頻域上的表達形式為該輸入信號的頻譜，頻譜反映了輸入信號的頻率構成，傅立葉變換是一種分析與處理一維信號的重要手段。當一個一維信號f?(x)滿足狄里赫萊條件，即f?(x)滿足：①具有有限個間斷點；②具有有限個極值點；③絕對可積＜，則此函數(shù)一定存在傅立葉變換對(傅立葉變換和逆變換)。在實際應用中，這些條件一般總是可以滿足的。一維傅立葉變換定義為逆變換定義為式中，j為虛數(shù)單位，稱為時(空)域變量，稱為頻域變量。(3.1)(3.2)一般情況下，函數(shù)f?(x)的傅立葉變換是一個復數(shù)量，可以表示為其中，為f?(x)的傅立葉幅度譜，為相位譜。E(u)?=?F2(u)為f?(x)的能量譜。(3.3)一維傅立葉變換很容易推廣到二維的情況，如果二維函數(shù)f?(x，y)滿足狄里赫萊條件，則它的二維傅立葉變換對為式中，x、y為時(空)域變量，u、v為頻域變量。(3.4)(3.5)類似于一維傅立葉變換，二維傅立葉變換的復數(shù)形式為

(3.6)其中，為f?(x,y)的傅立葉幅度譜，為相位譜，E(u，v)?=?|F(u，v)|2?=?R2(u，v)?+?I2(u，v)為能量譜。幅度譜表明了各正弦分量出現(xiàn)的多少，而相位譜信息表明了各正弦分量在圖像中出現(xiàn)的位置。對于整幅圖像來說，只要各正弦分量保持不變，幅值就顯得不是很重要。所以大多數(shù)時候使用濾波器都只影響幅值，而幾乎不改變相位信息。3.1.2離散傅立葉變換(DFT)連續(xù)傅立葉變換不適用于計算機的處理，而離散傅立葉變換由于其輸入和輸出都是離散值，故方便計算機的計算，且可以用快速傅立葉變換進行計算，以提高計算速度。本節(jié)介紹一維和二維的離散傅立葉變換。設{f?(0)，f?(1)，f?(2)，…，f?(N?-?1)}為一維信號f?(x)的N個采樣，由式(3.1)和式(3.2)可知，則可將離散傅立葉變換對定義為式中，x、u?=?0，1，2，…，N?-?1。式(3.8)中的系數(shù)1/N也可以放在式(3.7)中，有時也可在傅立葉正變換和逆變換前分別乘以1/，這是無關緊要的，只要正變換和逆變換前系數(shù)乘積等于1/N即可。(3.7)(3.8)式(3.7)和式(3.8)需要說明以下幾個概念：

(1)?f?(x)和F(u)都是離散序列。{f?(0)，f?(1)，f?(2)，…，f?(N?-?1)}表示取自相應連續(xù)函數(shù)的任意N個等間隔采樣值；當u?=?0，1，2，…，N?-?1時，F(xiàn)(u)?=?F(u0+uΔu)的值對應于在u0，u0+Δu，u0+2Δu，…，u0+(N-1)Δu處的連續(xù)變換的采樣值。

(2)頻域采樣間隔Δu與空域采樣間隔Δx的關系：Δu=1/(N·Δx)。

(3)?DFT總是存在的，不必考慮連續(xù)傅立葉變換的絕對可積條件。

(4)?DFT的f?(x)和F(u)都是周期性函數(shù)(周期為N)。在實際應用中，取一個周期，則f?(x)和F(u)是長度為N的有限序列。

根據(jù)式(3.7)，一維DFT可以表達為如下的矩陣形式：(3.9)式中，w=exp(-j2π/N)，稱為旋轉因子，wl×n表示w的l×n次方。常用一維DFT的核矩陣W是對稱矩陣，設G為一維信號矩陣，可將式(3.9)所示的一維離散傅立葉變換(DFT)用矩陣的形式表示為F?=?WG。同樣，很容易將一維離散傅立葉變換推廣到二維情況。二維離散傅立葉變換對定義為(3.10)(3.11)式中：x、u?=?0，1，2，…，M?-?1，y、v?=?0，1，2，…，N?-?1。與一維離散傅立葉變換類似，系數(shù)1/MN可以在正變換或逆變換中，也可以在正變換和逆變換前分別乘以系數(shù)1/，只要兩式系數(shù)的乘積等于1/MN即可。同樣，二維傅立葉變換也可以寫成矩陣形式：F?=?WgW

(3.12)其中，W為傅立葉變換核矩陣，g為二維信號矩陣。3.1.3二維傅立葉變換的性質(zhì)二維離散傅立葉變換具有以下一些重要的性質(zhì)(設函數(shù)f1(x，y)，f2(x，y)存在傅立葉變換)。

1.線性性質(zhì)設f1(x，y)，f2(x，y)的傅立葉變換分別為F1(u，v)，F(xiàn)2(u，v)，則函數(shù)af1(x，y)?+?bf2(x，y)的傅立葉變換為aF1(u，v)?+?bF2(u，v)，其中a、b為常數(shù)。

2.比例性質(zhì)設f(x，y)的傅立葉變換為F(u，v)，a、b為標量，則的傅立葉變換為aF(u，v)；f(ax，by)的傅立葉變換為。這說明空間比例尺展寬，對應于頻域比例尺壓縮。

3.可分離性由二維DFT的公式可知，正變換核可以分解成只含ux的和vy的的兩個指數(shù)函數(shù)的乘積，于是相應的二維DFT可以分離成連續(xù)兩次一維DFT(如圖3.1所示)：(3.13)(3.14)圖3.1用兩次一維DFT計算二維DFT

4.平移性質(zhì)設f?(x，y)的傅立葉變換為F(u，v)f?(x，y)F(u，v)，則根據(jù)傅立葉變換的定義可有平移性的表達：(3.15)(3.16)上式的物理意義：指數(shù)項乘f?(x，y)并取其乘積的變換，使頻率平面的原點移到點(u0，v0)；指數(shù)項乘F(u，v)并取其乘積的反變換，使空間平面的原點移到點(x0，y0)。同時應該注意到f?(x，y)的移動并不影響它的傅立葉變換的幅度。若取u0?=?v0?=?N/2，由于，從而，

(3.17)此時將f?(x，y)的傅立葉變換的原點移到相應N?×?N頻率方陣的中心。圖3.2給出這一過程的示意圖。圖3.2傅立葉頻譜平移示意圖

5.周期性傅立葉變換的正反變換x軸具有周期性，y軸具有N周期性，應用中只需取一個周期。在空間域中，f?(x，y)也有相似的性質(zhì)。

(3.18)(3.19)

6.共軛對稱性如果F(u，v)是f?(x，y)的傅立葉變換，F(xiàn)*(-u，-v)是F(u，v)變換的共軛函數(shù)，則有：F(u，v)?=?F*(-u，-v)(3.20)|F(u，v)|?=?|F(-u，-v)|(3.21)

7.旋轉不變性設f?(x，y)F(u，v)，且f?(x，y)和F(u，v)的極坐標形式分別為f?(r，θ)和F(ρ，φ)，則無論在連續(xù)的還是離散的傅立葉變換對中均有

(3.22)也就是說，如果f?(x，y)被旋轉θ°，則F(u，v)也被旋轉了同一角度；如果F(u，v)被旋轉θ°，則f?(x，y)也被旋轉了同一角度。圖3.3為這一過程的示意圖。圖3.3離散傅立葉變換的旋轉不變性

8.均值二維離散函數(shù)f?(x，y)的平均值定義為。將u?=?v?=?0代入式(3.10)，得從而

(3.23)

9.卷積定理兩個二維函數(shù)f?(x，y)和g(x，y)的卷積運算定義如下：

(3.24)其中，符號“”表示卷積運算。設f?(x，y)，g(x，y)的傅立葉變換分別為F(u，v)和G(u，v)，根據(jù)卷積和傅立葉變換的定義，傅立葉變換的卷積定理如下：

(3.25)(3.26)式(3.25)和式(3.26)表明時域(或空域)中的卷積等價于頻域的乘積。

10.相關定理兩個二維函數(shù)f?(x，y)和g(x，y)的相關函數(shù)定義如下：

(3.27)其中，符號“?！北硎鞠嚓P運算。則根據(jù)相關函數(shù)和傅立葉變換的定義可知傅立葉變換的另一個重要性質(zhì)——相關定理：··(3.28)(3.29)

11.微分性質(zhì)

f?(x，y)空間域上的微分f?(x，y)對應的傅立葉變換為(j2πu)m(j2πv)nF(u，v)。二元函數(shù)f?(x，y)的拉普拉斯(Laplacian)算子定義為

(3.30)按照二維傅立葉變換的定義可得拉普拉斯算子的傅立葉變換為

(3.31)其傅立葉變換等于傅立葉譜乘以u2?+?v2項，相當于傳遞函數(shù)隨頻率平方增加的線性系統(tǒng)。

12.?Parseval定理

Parseval定理也稱為能量保持定理。如果F(u，v)是f?(x，y)的傅立葉變換，那么有下式成立：

(3.32)這個性質(zhì)說明傅立葉變換前后的能量保持不變。3.1.4快速離散傅立葉變換離散傅立葉變換計算量非常大，運算時間長。可以證明其運算次數(shù)正比于離散點數(shù)N的平方，即N2，特別是當較大時，其運算時間將迅速增長。因此，研究離散傅立葉變換的快速算法(FastFourierTransform，F(xiàn)FT)是非常必要的。FFT算法可以將離散傅立葉變換的運算次數(shù)降為NlbN。其基本思路如下：

(1)將圖像進行轉置；

(2)按行對轉置后的圖像矩陣做一維FFT，將變換后的中間矩陣再轉置；

(3)對轉置后的矩陣做一維FFT，最后得到的就是二維FFT。

FFT算法只能處理大小為2的冪次的矩陣(其他大小的矩陣可以采用其他非基2的混合算法)。3.1.5MATLAB實現(xiàn)與舉例數(shù)字圖像本質(zhì)上為數(shù)組，所以MATLAB提供的對數(shù)組進行傅立葉變換操作的fft2、ifft2等函數(shù)同樣可以對圖像實現(xiàn)傅立葉變換。圖像處理中用得比較多的是fft2和ifft2函數(shù)，基本調(diào)用方法為

Y=fft2(X)

Y=fft2(X,m,n)

Y=ifft2(X)

Y=ifft2(X,m,n)其中X是輸入圖像矩陣，Y為返回圖像矩陣，X和Y的大小相同。fft2和ifft2函數(shù)的第二種調(diào)用方法，按照m和n指定的值先對圖像進行剪切或補零后再進行傅立葉變換，返回矩陣的大小為m×n。

以下的代碼首先調(diào)用fft2計算圖像的傅立葉變換，然后對此頻譜進行濾除處理，最后再重構出圖像：

RGB=imread('autumn.tif');

I=rgb2gray(RGB);

J=fft2(I); %圖像快速傅立葉變換

figure,imshow(I)J=fftshift(J); %頻譜平移figure,imshow(log(abs(J)),[])J(abs(J)<5000)=0; %小信息濾除figure,imshow(log(abs(J)+eps),[])J=ifftshift(J);K=ifft2(J); %圖像傅立葉逆變換figure,imshow(K,[0255])運行結果如圖3.4所示。圖3.4圖像傅立葉變換3.2離散余弦變換(DCT)

3.2.1一維離散余弦變換余弦變換為傅立葉變換的特例。設一維離散序列f?(x)，x?=?1，2，…，N?-?1，以-1/2為中心反折，形成-N到-1的序列，與原序列合并形成長度為2N的偶序列，此時傅立葉變換的核函數(shù)exp(-j2πux/N)改變?yōu)閑xp(-j2π(x?+?1/2)u/2N)，按傅立葉變換的性質(zhì)，虛部為零不進行運算，核函數(shù)為余弦函數(shù)，等價于cos[(2x?+?1)uπ/2N]。根據(jù)傅立葉變換的定義，余弦正變換為(3.33)為保證每行正交向量模為1，對上式進行歸一化處理：

u?=?0，1，2，…，N?-?1(3.34)其中

(3.35)

k?=?0，1，2，…，N?-?1(3.36)一維余弦變換用矩陣可表示為FC?=?Cf

(3.37)其中核矩陣C為

(3.38)一維DCT的核矩陣是正交矩陣，但不是對稱矩陣。根據(jù)正交性，逆變換矩陣為(3.39)因此可見一維DCT的逆變換核與正變換核是相同的。因為余弦變換是傅立葉變換的特例，傅立葉反變換的核矩陣是W陣的共軛矩陣，對于余弦變換，共軛矩陣等于自身，因此(注意：除了行、列序號互換外，形式上與正變換完全一樣)

n?=?0，1，2，…，N?-?1(3.40)其矩陣形式為f?=?CTFC(3.41)3.2.2二維離散余弦變換很容易將一維DCT的定義推廣到二維DCT(思考：如何形成二維偶函數(shù)？先水平作對折鏡像，然后再垂直作對折鏡像)。數(shù)字圖像f?(m，n)可看成是一個M×N的矩陣，通過二維DCT，可以將圖像從空間域變換到DCT域。二維DCT正變換核為

(3.42)式中，C(u)和C(v)的定義同式(3.42)；x、u?=?0，1，2，…，M?-?1；y、v?=?0，1，2，…，N?-?1。設f(x，y)為M?×?N的數(shù)字圖像矩陣，則二維DCT定義如下：

(3.43)二維DCT逆變換定義如下：

(3.44)類似于一維矩陣形式的DCT，可以寫出二維DCT的矩陣形式：F?=?CfCT(3.45)二維DCT的逆變換的矩陣形式為f?=?CTFCC(3.46)余弦變換具有如下性質(zhì)：

(1)余弦變換為正交變換，即C?=?C*，C-1?=?CT；

(2)離散序列的余弦變換是DFT的對稱擴展形式；

(3)和傅立葉變換類似，余弦變換也存在快速變換；

(4)和傅立葉變換類似，余弦變換也具有將高度相關數(shù)據(jù)能量集中的優(yōu)勢。數(shù)字圖像處理中最為常用的DFT快速算法，需要進行算術運算，可將整幅圖像的信息用若干個系數(shù)很好地表達出來。DCT的快速算法，只要求實數(shù)運算，在高相關性圖像的處理中，它最接近最佳的KL變換，用于實現(xiàn)編碼和維納濾波。3.2.3MATLAB實現(xiàn)函數(shù)在MATLAB中，dctmtx函數(shù)用來創(chuàng)建離散余弦變換的核矩陣，調(diào)用方法如下：D=dctmtx(n)其中輸入變量n為double類型的整數(shù)，返回矩陣D也為double類型的。調(diào)用dctmtx函數(shù)創(chuàng)建n?×?n離散余弦變換核矩陣，利用該核矩陣可以直接計算余弦變換或者逆余弦變換。列向量A的離散余弦變換為D*A。A的逆余弦變換為D'A。如果A為方陣，則矩陣A的二維離散余弦變換可以計算為DAD'，這種計算方法比直接調(diào)用dct2函數(shù)速度要快，特別是需要計算大量低階矩陣的DCT時。這是因為這種計算方法只用計算一次D。根據(jù)式(3.38)，當N?=?4時，余弦變換核矩陣C為(3.47)根據(jù)式(3.38)，還可以計算出向量f?=?[1331]'的余弦變換：(3.48)根據(jù)式(3.45)，可求出圖像f1的余弦變換F1為(3.49)上述計算過程用MATLAB相關函數(shù)可以表達成

C=dctmtx(4)

f=[1331]';

F=C*f

f1=[1331;1331;1331;1331];

F1=C*f1*C'計算結果為C=0.5000 ?0.5000 ?0.5000 ?0.50000.6533 ?0.2706? -0.2706 -0.65330.5000 -0.5000 -0.5000 ??0.50000.2706 -0.6533 ?0.6533 -0.2706F=4.00000.0000-2.0000-0.0000F1=?8.0000?0.0000-4.0000-0.0000?0.0000?0.0000-0.0000-0.0000-0.0000-0.0000?0.0000?0.0000-0.0000-0.0000?0.0000?0.0000

MATLAB提供的對數(shù)組進行離散余弦變換操作的dct2、idct2等函數(shù)同樣可以對圖像實現(xiàn)離散余弦變換。圖像處理中用得較多的是dct2和idct2函數(shù)，基本調(diào)用方法分別如下：

B=dct2(A)

B=dct2(A,m,n)

B=dct2(A,[mn])

B=idct2(A)

B=idct2(A,m,n)

B=idct2(A,[mn])其中，A是輸入圖像矩陣，B為返回圖像矩陣，A和B的大小相同。fft2和ifft2函數(shù)的第二種調(diào)用方法是按照m和n指定的值對圖像進行剪切或者補零后進行余弦變換，返回矩陣的大小為m×n。第三種調(diào)用方法與第二種相同。

例如，下面的代碼首先調(diào)用dct2計算圖像的余弦變換，然后對此變換結果進行濾除處理，最后再重構出圖像：

RGB=imread('autumn.tif');

I=rgb2gray(RGB);

J=dct2(I);

figure,imshow(I)

figure,imshow(log(abs(J)),[])J(abs(J)<10)=0;K=idct2(J);figure,imshow(log(abs(J)),[])figure,imshow(K,[0255])運行結果如圖3.5所示。圖3.5圖像余弦變換3.3射線投影積分變換

沿著某條射線對圖像進行積分得到的投影變換為投影積分變換，也稱Radon變換。從此變換矩陣恢復出原始圖像稱為斷層影像重建。這種變換方法在醫(yī)學領域應用廣泛，如CT。根據(jù)投影線之間的關系可以將此類投影變換分為兩類：平行線投影(圖3.6)和扇形束投影變換(圖3.7)。下面分別介紹這兩類投影的正逆變換，并介紹MATLAB實現(xiàn)此變換的函數(shù)和示例。圖3.6平行投影

圖3.7扇形束投影3.3.1平行線投影變換和逆變換

1.定義如圖3.6所示，平行線投影變換可定義為圖像函數(shù)沿著直線的線積分。二維函數(shù)f?(x，y)的投影就是一系列線積分。而圖像重建過程即將圖像函數(shù)從許多不同角度下沿直線的線積分所產(chǎn)生的投影數(shù)據(jù)進行逆Radon變換，從而得到截面函數(shù)。MATLAB中的Radon函數(shù)可計算不同角度下的線積分。根據(jù)定義，圖3.8所示的矩形圖像f?(x，y)沿著豎直方向的直線的積分就是函數(shù)f?(x，y)到x軸的投影，它沿著水平方向的直線的積分就是函數(shù)f?(x，y)到y(tǒng)軸的投影。圖3.8簡單圖像沿著豎直和水平方向的投影如圖3.9所示，f?(x，y)的Radon變換就是沿著平行于軸的直線的線積分：

(3.50)其中

(3.51)圖3.9Radon變換沿任意方向投影計算

Radon逆變換也就是根據(jù)Radon變換的投影數(shù)據(jù)重建原始圖像，常用于斷層影像重建中。在x-ray斷層影像重建中，變換數(shù)據(jù)通過指定物體的不同角度射線衰減得到。原始圖像可認為是樣本的某個截面；變換數(shù)據(jù)利用特定用途的硬件設備收集得到，利用這些數(shù)據(jù)，調(diào)用iradon函數(shù)可以重構出原始圖像。這種方法可以在不破壞活體或非透明樣本的情況下獲取其內(nèi)部圖像。在平行束幾何中，每個投影數(shù)據(jù)是通過在某一角度下的一系列直線上積分得到的。平行束幾何在x射線吸收層析照相術中的應用示意圖如圖3.10所示。需要注意的是這里的射線發(fā)射器數(shù)和接受器數(shù)是相等的。每個接受傳感器測量其對應的射線發(fā)射器發(fā)送的射線，射線的衰減給出了待測物體的綜合密度、質(zhì)量的某種度量。這個測量值對應著Radon變換中的線積分。圖3.10Radon逆變換計算變換成像法中傅立葉切片定理聯(lián)系了Radon變換與傅立葉變換之間關系。該定理將投影數(shù)據(jù)與物體截面的二維傅立葉變換關聯(lián)起來。傅立葉切片定理的含義是將平行投影的一維傅立葉變換等同于原始物體的二維傅立葉變換的一個切片，如圖3.11所示。圖3.11傅立葉切片定義圖像重建過程實際上就是逆Radon變換過程，具體步驟如下：

(1)投影的傅立葉變換為頻域響應某個傾斜角度的一維濾波器濾波；

(2)將濾波后的結果進行傅立葉逆變換；

(3)最后進行反投影得到重建的圖像。濾波運算可以在空域或頻域進行，從而有不同的重建算法。

2.?MATLAB實現(xiàn)

MATLAB提供radon函數(shù)來計算圖像在指定角度上的Radon變換，其基本調(diào)用格式如下：[R,xp]=radon(I,theta)其中，I為待變換的圖像，R的各行為返回角度參數(shù)theta中各方向上的radon變換值；xp表示沿x'?軸相應的坐標值，I的中心定義為floor((size(I)?+?1)/2)。這個點位于x'?=?0的x'軸上。例如，下面的代碼計算并繪制單一正方形沿著0°和45°的Radon變換。xp對任意角度的變換都保持不變。

I=zeros(100,100);

I(25:75,25:75)=1;

imshow(I)

[R,xp]=radon(I,[045]);

figure;plot(xp,R(:,1));title('R_{0^o}(x\prime)')

figure;plot(xp,R(:,2));title('R_{45^o}(x\prime)')運行結果如圖3.12所示。圖3.12Radon變換例程一結果圖像沿著一系列角度的Radon變換通常一起顯示以反映圖像的某種屬性。例如，以下的代碼計算上述的矩形沿著0°到180°的，每增加1°的Radon變換值。

theta=0:180;

[R,xp]=radon(I,theta);

imagesc(theta,xp,R);

title('R_{\theta}(X\prime)');

xlabel('\theta(degrees)');

ylabel('X\prime');

set(gca,'XTick',0:20:180);

colormap(hot);

colorbar運行結果如圖3.13所示。圖3.13Radon變換例程二結果

MATLAB提供iradon函數(shù)計算Radon逆變換，也即從Radon變換R中重構出圖像IR，其基本調(diào)用格式為IR=iradon(R,theta)

iradon函數(shù)使用濾波變換(filterbackprojection)算法來計算Radon逆變換。這個算法基于投影數(shù)據(jù)R中的列重構出原始圖像I的近似值。更準確的重構結果可以使用更多的投影數(shù)據(jù)得到。隨著投影數(shù)據(jù)的增多(theta的個數(shù))，重構圖像IR可以更準確地描述原始圖像I。theta向量必須包含單調(diào)增的角度值，且角度增量Dtheta必須相等。當角度增量Dtheta已知，可以將此量作為iradon函數(shù)輸入變量theta的替代，也就是IR=iradon(R,Dtheta)。濾波變換算法對投影數(shù)據(jù)R進行濾波，然后再重構出圖像。在某些情況下，投影數(shù)據(jù)中存在噪聲。為了消除高頻噪聲，可以利用某種窗函數(shù)對投影數(shù)據(jù)進行處理。大多數(shù)窗函數(shù)濾波器可以在iradon函數(shù)中使用。更多的說明參見iradon函數(shù)。

iradon函數(shù)允許用戶指定歸一化的頻率D，大于此頻率時，濾波器的輸出為零。D的取值為[0，1]。例如，以下代碼首先調(diào)用Radon函數(shù)計算某一圖像的Radon變換，然后調(diào)用iradon函數(shù)重構出此圖像：P=phantom(128);theta=0:179;[R,xp]=radon(P,theta);IR=iradon(R,theta,'nearest','Hann');imshow(P)figure,imagesc(theta,xp,R);title('R_{\theta}(X\prime)');xlabel('\theta(degrees)');ylabel('X\prime');set(gca,'XTick',0:20:180);colormap(hot);colorbarfigure,imshow(IR)運行結果如圖3.14所示。圖3.14Radon逆變換例程結果

3.應用實例例3.1直線檢測。

Radon變換與計算機圖像操作中的Hough變換類似?？梢哉{(diào)用radon函數(shù)實現(xiàn)某一種Hough變換以檢測直線。本例利用Radon變換檢測直線，具體實現(xiàn)步驟如下：

(1)利用edge函數(shù)計算圖像的邊緣以將圖像轉化成二值圖像。

I=fitsread('solarspectra.fts');

I=mat2gray(I);

BW=edge(I);

imshow(I),figure,imshow(BW)運行結果如圖3.15所示。圖3.15待檢測圖像(2)計算獲得二值圖像的Radon變換。

theta=0:179;

[R,xp]=radon(BW,theta);

figure,imagesc(theta,xp,R);colormap(hot);

xlabel('\theta(degrees)');ylabel('x\prime');

title('R_{\theta}(x\prime)');

colorbar運行結果如圖3.16所示。圖3.16待檢測二值圖像的Radon變換

(3)尋找Radon變換矩陣中的峰值點，這些峰值點對應著原始圖像中的直線。圖3.17顯示，Radon變換矩陣R中最強烈的點的坐標為θ?=?1°和x'?=?-80。圖3.17中虛線方框內(nèi)的線段是垂直于這個方向且過x'?=?-80的直線。Radon變換幾何以黑色顯示。與虛線框內(nèi)的線平行的直線在Radon變換圖像中也出現(xiàn)在θ?=?1°的位置。同樣，垂直于這個直線的直線在Radon變換矩陣中的位置為θ?=?91°。圖3.17Radon變換峰值位置示意圖例3.2頭顱圖像重構。圖3.18原始圖像

(1)讀入原始圖像，原始圖像如圖3.18所示。

P=phantom(256);

imshow(P)

圖3.18原始圖像

(2)計算圖像的Radon變換。本例程計算三種不同角度采樣間隔下的Radon變換，以比較分析結果。

theta1=0:10:170;%角度參數(shù)間隔為10°

[R1,xp]=radon(P,theta1);

num_angles_R1=size(R1,2);

theta2=0:5:175;%角度參數(shù)間隔為5°

[R2,xp]=radon(P,theta2);

num_angles_R2=size(R2,2)

theta3=0:2:178;%角度參數(shù)間隔為2°

[R3,xp]=radon(P,theta3);

num_angles_R3=size(R3,2)這里以角度參數(shù)間隔2°的情形為例，繪制出圖像的Radon變換結果，如圖3.19所示。

figure,imagesc(theta3,xp,R3)

colormap(hot)

colorbar

xlabel('ParallelRotationAngle-\theta(degrees)');

ylabel('ParallelSensorPosition-x\prime(pixels)');圖3.19Radon變換結果

(3)利用iradon重構圖像。分別計算上述三個角度參數(shù)下的重構圖像，代碼如下：

output_size=max(size(P));

dtheta1=theta1(2)-theta1(1);%角度參數(shù)10°時的重建

I1=iradon(R1,dtheta1,output_size);

figure,imshow(I1)

dtheta2=theta2(2)-theta2(1);%角度參數(shù)5°時的重建

I2=iradon(R2,dtheta2,output_size);

figure,imshow(I2)

dtheta3=theta3(2)-theta3(1);%角度參數(shù)2°時的重建

I3=iradon(R3,dtheta3,output_size);

figure,imshow(I3)圖3.20為運行結果圖。由圖可看出，隨著采樣間隔的減小，重構結果越來越精細。圖3.20Radon逆變換重構結果3.3.2扇束變換和逆變換

1.定義和實現(xiàn)扇束變換計算某一圖像沿著特定方向的投影數(shù)據(jù)，如圖3.21所示。二維圖像f?(x，y)的投影就是沿著一系列直線的積分。MATLAB用fanbeam函數(shù)計算沿著從單點出發(fā)的直線的線積分，這些投影直線匯交于一點，形成了扇形。為了能夠表達一個圖像，必須使源繞著圖像中心旋轉以形成不同角度的積分值。

如圖3.21所示，扇形投影幾何的傳感器有兩種排列方式，分別為弧形排列和直線形排列。源中心點到圖像中心的距離D、弧形排列中相連傳感器之間的夾角、線形排列中傳感器相連之間的距離皆可控。圖3.21fanbeam計算幾何示意圖

MATLAB調(diào)用fanbeam函數(shù)計算扇束變換投影數(shù)據(jù)，其基本調(diào)用格式為F=fanbeam(I,D)以待計算圖像I和源中心與圖像旋轉中心之間的距離D為輸入變量，fanbeam函數(shù)利用圖像的大小和fanbeam參數(shù)來決定束數(shù)。fanbeam函數(shù)的可選參數(shù)有FansensorGeometry、FansensorSpacing和FanRotationIncrement。

FansensorGeometry選項指定接受傳感器的排列方式，可為“arc”或者“l(fā)ine”，其中“arc”為缺省值。當為“arc”時，fanbeam