卷积反投影重建(2维).docxVIP

卷积反投影图像重建1反投影重建基本介绍设待重建图像为，它的二维傅氏变换为。根据中心切片定理，可通过在不同视角下的投影的一维傅氏变换求得。即：待建图像：(1.1)因为，所以有：同时： (1.2)先来看该式的第二个积分：(1.3)式中：(1.4)式(3.10)的物理意义是投影经过传递函数为的滤波器后得到的修正后的投影在满足时的值。将(3.11)代入(3.8)，得到： (1.5)称为滤波反投影方程，其物理意义是经过给定点的所有滤波后的投影在~范围内的累加—反投影重建，得出点的像素值。可见，滤波（卷积）反投影算法的具体包含三大步：(1) 把在固定视角下测得的投影经过滤波，得到滤波后的投影；(2) 对每一个，把反投影于满足的射线上的所有各点；(3) 将步骤(2)中的反投影值对所有进行累加（积分），得到重建后的图像。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 scansrand(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(:)));endR1 = log(I0 ./ max(yi,0.01)); % noisy sinogramR1=max(R1,0); % 显示ff=2;uu=22000;v=ff*exp(R1/uu);subplot(2,2,2)imagesc(R1);title(256*360有噪声平行投影);colormap(gray)colorbarQ=reshape(R1,256,360);效果图如图-2: 图-22.5 滤波器的选择与设计最基本的从投影重建图像的滤波器：1971年提出的R-L重建滤波器(下图中实线表示)。图-3 R-L 滤波器示意图空域表达式为:(2.1)其中B为截至频率。若，d为空间采样间隔，可解出离散的滤波器空间脉冲响应： (2.2)MATLAB 代码如下:g=-(Nd/2-1):(Nd/2);for i=1:256 if g(i)==0 hl(i)=1/(4*d^2); else if mod(g(i),2)==0 hl(i)=0; else hl(i)=(-1)/(pi^2*d^2*(g(i)^2)); end endendk=Nd/2:(3*Nd/2-1);2.6 重建图像代码设计:for m=1:Nt pm=Q(:,m); u=conv(hl,pm); pm=u(k); Cm=((