卷积反投影重建(二维)
![卷积反投影重建(二维)](https://img.360docs.net/imga1/04us18up20tzladx1e8x-11.webp)
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
卷积反投影图像重建
1 反投影重建基本介绍
设待重建图像为),(y x a ,它的二维傅氏变换为^
1
2
(,)(,)A A ωωρφ=。根据中心切片定理,^
(,)A ρφ可通过),(y x a 在不同视角φ下的投影()r p x φ的
一维傅氏变换求得。即:
待建图像:
12^
1212()1212
2
^2cos()0
2cos()02cos()0(,)(,)(,)1(,)4(,)(,)(,)i x y i r i r i r a r a x y F A A e d d A e d d P e d d d P e d ωωππρθφπ
πρθφπ
πρθφθωωωωωωπρφρρφ
ρφρρφ
φρρφρ
-∞∞
+-∞-∞
∞--∞
∞
--∞
∞
--∞======⎰⎰
⎰⎰⎰
⎰
⎰⎰
[]
(1.1)
因为cos()r x r θφ=-,所以有:
122(cos sin )22cos()r x y x y x r ωωπρφφπρπρθφ+=+==- 同时:
12d d J d d ωωρφ=
11222//2cos 2sin 4//2sin 2cos J ωρωφπφ
πρφ
πρωρωφπφ
πρφ
∂∂∂∂-=
=
=∂∂∂∂ (1.2)
先来看该式的第二个积分:
[]
22cos()
cos()cos()cos()(,)(,)|()(,)|(,)|cos(),r r r r i x i r x r r r x r r x r P e
d P
e d h x p x g x g r πρπρθφθφθφθφρρφρρρφρ
φφθφφ∞
∞
-=--∞
-∞
=-=-==*==-⎰
⎰ (1.3)
式中:(,)()(,)r r r g x h x p x φφ=*
(1.4)
^
121(,)(,)()()(,)r
A A F p x P
P φφωωρφρρφ⎡⎤====⎣⎦
式(3.10)的物理意义是投影(,)r p x φ经过传递函数为1[()]r F h x ρ=的滤波器后得到的修正后的投影(,)r g x φ在满足cos()r x r θφ=-时的值。将(3.11)代入(3.8),得到:
^
0(,)[cos(),]a r g r d π
θθφφφ=-⎰ (1.5)
称为滤波反投影方程,其物理意义是经过给定点(,)r θ的所有滤波后的投影在0φ=~π范围内的累加—反投影重建,得出(,)r θ点的像素值。
可见,滤波(卷积)反投影算法的具体包含三大步:
(1) 把在固定视角i φ下测得的投影(,)r p x φ经过滤波,得到滤波后的投影(,)r g x φ;
(2) 对每一个i φ,把(,)r i g x φ反投影于满足cos()r i x r θφ=-的射线上的所有各点(,)r θ;
(3) 将步骤(2)中的反投影值对所有0<φπ≤进行累加(积分),
得到重建后的图像。 2 重建流程
2.1 首先我们利用phantom ()函数产生一个头部幻影图像,用以检测二维重建算法,代码如下:
I=phantom(256); subplot(2,2,1) imshow(I,[]);
title('256*256原始图像'); 效果图如图1所示,
图1
为一个大椭圆和几个小椭圆。
2.2 初始参数设置
重建采用的是平移加旋转的扫描方式,射线源在某一角度下水平移动,将物体全部照射后旋转一角度,如此重复,在这个过程中探测器相应地运动以接收X 射线。根据此原理,将重建程序的初始参数设置如下
:
[N,N]=size(I);
z=2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3;% radon变换默认平移点数/角度
Nt=360; % 角度采样点数
Nd=N; % 平移数
x=pi/180; % 角度增量
d=N/Nd; % 平移步长
theta = 1:Nt;
a=zeros(N);
2.3 产生无噪声投影数据
[R,xp] = radon(I,theta);
e=floor((z-Nd)/2)+2;
R=R(e:(Nd+e-1),:);
R1=reshape(R,256,360);
radon(I,theta)产生I投影,默认z点/角度,即使指定N点也是z 点.所以为避免重建图像放大或缩小,下面计算取投影时需补偿,补偿量e,如对256的图像,补偿为55,即pm的第55个点作为计算用的第一个投影。
2.4 添加噪声并将所有噪声平行投影进行显示
[mm,nn]=size(R1);
di=lognrnd(0,0.15,mm,nn);
R1= 10*(R1-min(R1(:)))/( max(R1(:))-min(R1(:)));
I0 = 1.5e5; % incident photons; decrease this for simulating "low dose" scans
rand('state', 0), randn('state', 0);
yi= poissrnd(I0 * di.*exp(-R1))+3*randn(size(R1));
if any(yi(:) == 0)
warn('%d of %d values are 0 in sinogram!', ...
sum(yi(:)==0), length(yi(:)));
end
R1 = log(I0 ./ max(yi,0.01)); % noisy sinogram
R1=max(R1,0);
% 显示
ff=2;
uu=22000;
v=ff*exp(R1/uu);
subplot(2,2,2)
imagesc(R1);
title('256*360有噪声平行投影');
colormap(gray)
colorbar
Q=reshape(R1,256,360);
效果图如图-2: 图-2 2.5 滤波器的选择与设计