1.算法描述
從投影重建物體的截面圖像是圖像處理中非常重要的技術(shù)此技術(shù)在物體的無損傷性檢測(cè)其內(nèi)部缺陷的應(yīng)用中能起很大作用從投影重建圖像的技術(shù)早在20世紀(jì)中期就已經(jīng)制成常規(guī)醫(yī)療診斷設(shè)備的商品1917年奧地利數(shù)學(xué)家J.Radon發(fā)表的論文證明了二維物體或三維物體可以從許多投影來重建其內(nèi)部數(shù)據(jù)英國的EMI公司中央研究所1972年由G.N.Hounsfield設(shè)計(jì)出X射線斷層分析儀此后于1974年5月在加拿大蒙特利爾第一次國際會(huì)議上正式命名為Computerized Tomography簡稱XCT或CT當(dāng)時(shí)這臺(tái)儀器可以獲得人體的各個(gè)部位鮮明清晰的斷層圖像另外美國塔夫茨大學(xué)A.M.Cormack在1963年就提出了精確的非迭代的級(jí)數(shù)展開法用于圖像重建的算法中并指出在診斷醫(yī)學(xué)中圖像重建的可能性以上兩人由于對(duì)CT研制工作的開創(chuàng)性貢獻(xiàn)獲得了1979年諾貝爾醫(yī)學(xué)獎(jiǎng)缺菌。
計(jì)算機(jī)層析成像技術(shù)（Computed Tomography，簡稱CT）是利用具有一定能量的射線源（X射線核芽，γ射線）對(duì)物體進(jìn)行斷層掃描亡脑，并根據(jù)物體外部的探測(cè)器獲得的物理量（指物質(zhì)對(duì)射線的衰減系數(shù)）生成的一維投影數(shù)據(jù)诅需，通過特定的重建算法渔呵，得到所掃描斷層的二維圖像责掏。CT斷層圖像具有無影像重疊、空間和密度分辨率高轩拨、可直接進(jìn)行數(shù)字化處理等優(yōu)點(diǎn)践瓷，通過近十幾年的發(fā)展已成為非接觸無損檢測(cè)的主流技術(shù)，是關(guān)鍵部件檢測(cè)亡蓉、機(jī)械仿型設(shè)計(jì)晕翠、安全檢查等方面強(qiáng)有力的手段，并廣泛應(yīng)用于航空砍濒、航天淋肾、機(jī)械、汽車梯影、船舶巫员、公安等領(lǐng)域庶香。此外甲棍，在醫(yī)學(xué)領(lǐng)域，醫(yī)用CT機(jī)已成為疾病診斷的重要輔助工具赶掖，對(duì)病灶的位置及病變程度的良好再現(xiàn)已使其成為醫(yī)學(xué)影像領(lǐng)域里的主要醫(yī)療器械感猛。CT機(jī)有三個(gè)組成部分。一是數(shù)據(jù)采集系統(tǒng)奢赂，二是圖像重建（主要是圖像重建算法）陪白，三是圖像后處理及顯示系統(tǒng)。CT圖像重建算法主要分為兩類膳灶，即變換法和級(jí)數(shù)展開法咱士。變換法包括濾波反投影算法立由、傅立葉變換法以及ρ濾波法（rho-filtered layergrams），小波變換法序厉；級(jí)數(shù)展開法主要是代數(shù)重建算法锐膜。濾波反投影是比較常用的變換算法，它具有速度快弛房，空間和密度分辨率好的優(yōu)點(diǎn)道盏。
濾波反投影法是目前CT圖像重建領(lǐng)域使用得最為廣泛的一種算法，其有效性和準(zhǔn)確性也在臨床上得到了驗(yàn)證文捶。
通常在計(jì)算機(jī)中荷逞，采用的離散反投影算法，因此離散情形下的卷積反投影算法可以如下描述：
其中O(x粹排，y)代表原始圖像种远，R(x，y)代表重建圖像恨搓。當(dāng)然孤立地去看相對(duì)誤差并不是一個(gè)好的衡量標(biāo)準(zhǔn)院促，它必須同人們的視覺評(píng)估結(jié)合起來進(jìn)行考慮。
shepp-Logan模型:
用S-L頭模型進(jìn)行計(jì)算機(jī)仿真研究的主要好處之一是可以獲得該模型投影數(shù)據(jù)的解析表達(dá)式斧抱。一個(gè)中心位置在原點(diǎn)且未經(jīng)旋轉(zhuǎn)的橢圓常拓，其長軸與X軸重合，短軸與Y軸重合辉浦。假設(shè)橢圓內(nèi)的密度為r弄抬、橢圓外密度為零，則該橢圓圖像可用以下方程表示：
式中A宪郊，B分別為橢圓長軸與短軸的長度掂恕。
在研究從投影重建圖像的算法時(shí)，為了比較客觀的評(píng)價(jià)各種重建算法的有效性弛槐，人們常選用公認(rèn)的Sheep-Logan 頭模型作為研究對(duì)象懊亡。該模型由10個(gè)位置、大小乎串、方向店枣、密度各異的橢圓組成，象征一個(gè)腦斷層圖像叹誉。
該模型的主要參數(shù)如下表所示：
多尺度全局重建算法計(jì)算流程:
2.仿真效果預(yù)覽
matlab2022a仿真結(jié)果如下：
3.MATLAB核心程序
function filtPR=lvbo(PR);
n = size(PR,1);
sideSize = n;
a = 1;
[Length, Count] = size(PR);
w = [-pi:(2*pi)/Length:pi-(2*pi)/Length];
rn1 = abs(2/a*sin(a.*w./2));
rn2 = sin(a.*w./2);
rd = (a*w)./2;
r = rn1*(rn2/rd)^2;
f = fftshift(r);
for i = 1:Count
IMG = fft(PR(:,i));
fimg = IMG.*f';
g(:,i) = ifft(fimg);%利用FFT,IFFT進(jìn)行濾波
end
filtPR= real(g);
............................................................
I=phantom(256); ?%生產(chǎn)頭部模型圖
figure(1);
imshow(I); ??????%顯示圖像
IMG=double(I); ??%雙精度顯示
[cod_a,cod_h,cod_v,cod_d,map]=wtest(IMG);%一層小波變換
figure(2);
subplot(221)%顯示分解系數(shù)1
imshow(cod_a,map);
subplot(222)%顯示分解系數(shù)2
imshow(cod_h,map);
subplot(223)%顯示分解系數(shù)3
imshow(cod_v,map);
subplot(224)%顯示分解系數(shù)4
imshow(cod_d,map);
Z1=idwt2(cod_a,cod_h,cod_v,cod_d,'db1');%重建
figure(3);
imshow(Z1,map);
for i=1:256
for j=1:256
error(i,j)=(I(i,j)-Z1(i,j))^2/I(i,j)^2;%誤差求解
end
end
for i=1:256
for j=1:256
if I(i,j)==0
error(i,j)=255; ???
end
end
end
for i=1:256
for j=1:256
if error(i,j)>255;
error(i,j)=255; ???
end
end
end
