平行束滤波反投影1100500121 赵伟伦 准备知识:一维Fourier 变换：dt et f f f F t i ⎰+∞∞--⋅==πωω2)()(~)( 一维逆Fourier 变换： ωωπωd e f f F x f t i ⎰+∞∞--⋅==21)(~)~()( 且有：)~(~),(11f F F f f F F f --⋅=⋅=重要的性质：（卷积特性）)(~)(~)*(ωωgf g f F ⋅=； )(~)(~)(ωωgf g f F *=⋅ 二维Fourier 变换： dX e x x f f f F x x i R ),(),(22121221212),(),(~)(⋅-⎰==ωωπωω； 逆二维Fourier 变换： Ω==⋅-⎰d e f f F x x f x x i R ),(),(221122121212),(~)~(),(ωωπωω； 中心切片定理：),)(ˆ()(2ϕωωfF f F r =Φ, 其中),(ˆϕr f 是),(21x x f 的Radon 变换： 解释：一个二元函数的Radon 变换关于r 的一维Fourier 变换与这个二元函数的二维Fourier 变换形式相等。

滤波反投影：思路：)(),(121f F F x x f ⋅=-()()[][]ϕϕωωϕωϕωϕωωϕωϕωϕωωωϕωωϕωϕωωϕωϕωωωϕωωωππωωππωωππωωππωωπd r f F r d fF F d d e fF x x r d d e fF d d e f F d d e f d d e f F X r x x r r r r i r x x i r x x i rx x i x x i R Φ⋅=-Φ⋅=-∞+∞-⋅∞+∞-⋅∞+⋅∞+⋅*⇔=⋅⇔⇔Φ⋅=Φ=⇔⇔⇔⇔⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰)(H ),(ˆfourier fourier ),()(H ),)(ˆ(]),)(ˆ([),),(),(),(),)(ˆ(),)(ˆ()(~)(1),(1202121),(),(20),(),(2200),(),(2200221),(),(222121212121212121212变化变化等于函数点乘后的个函数的卷积的并根据卷积的性质：两设旋转角为为坐标映射到探测器上，设为用极坐标方式表示出来（把，可知），（由于中心切片定理)(),(~),(r H r f r G *=ϕϕ)(r H 是滤波器总结：ϕϕϕωωϕωππωπd r H r fd def F X f X r X r r i r Φ⋅=Φ⋅=+∞∞-⎰⎰⎰=⎥⎦⎤⎢⎣⎡⋅=)(*),(ˆ),)(ˆ()(020 解释为：投影数据),(ˆϕr f 先进行滤波)(*),(ˆr H r f ϕ 在对滤波数据进行投影ϕϕπd r H r f X r Φ⋅=⎰)(*),(ˆ0简单例子：（大圆与小圆）通过已得到的正投影‘round.dat’经过滤波后，反投影后的图像。

正投影数据：滤波图像：反投影后的图像：总的滤波反投影过程：1，得到图像的正投影2，滤波（投影与滤波器卷积）3，反投影（一个像素所有角度的滤波函数值的和再乘以pi/views ）views r H r P d r H r P x x f viewsπϕϕϕϕπ∆=*=∑⎰=)(*),()(),(),(0021例子：GUI界面实现：修改前：程序讲述：主程序：读入一张512*512的黑白图片，展示图片；求正投，展示正投结果；求滤波后的结果并展示；求反投影后的的结果并展示。

view=360;bins=512;U=1.5; %U=view_radiusP=zeros(view,bins);clc;clear;A = imread('angle.jpg');width = 512;height = 512;A = im2double(A);figure,imshow(A)t1 = -1;t2 = 1;n= 2048;dt = (t2-t1)/n;bins = 512;view = 360;du = 3/(bins-1);P = zeros(view,bins);time = clock;for i = 1:viewphi = (i - 1)*(pi/180);for j = 1:bins%将视野范围512等分（定义512个接受器）r = (j - 1)*du - 1.5 + 0.5*du;%设定第一个接受器的与探测器中心的距离int = integral(r,phi,t1,t2,dt,width,height,A);P(i,j) = int;endendouttime = etime(clock,time);s_png2=sprintf ('angle.dat') ;s_png2=strcat ('D:\文件\ct\zwl\',s_png2) ;fp = fopen(s_png2,'wb');fwrite(fp,P','float64');fclose(fp);s_png = sprintf('angle.dat');%angles_png = strcat('D:\文件\ct\zwl\',s_png);fp = fopen(s_png,'r');P = fread(fp,view*bins,'float64');fclose(fp);P = reshape(P,bins,view);P = P';U=1.5;cellsize=2*U/bins;theta0=2*pi/view;figure, imshow(P);integ=RL_filter(P,view,bins,cellsize);figure, imshow(integ,[]);Image=back_projection_circle(integ,cellsize,view,theta0,U);figure,imshow(Image,[])s_png2=sprintf ('R_L_projection1.dat') ;s_png2=strcat ('',s_png2) ;fp = fopen(s_png2,'wb');fwrite(fp,Image','float64');fclose(fp);子函数1：备注：将图像积分区域转化为二维图像中的点的方法：function int = integral(r,phi,t1,t2,dt,width,height,A)int = 0;for t = t1:dt:t2x = r*cos(phi) - t*sin(phi);y = r*sin(phi) + t*cos(phi);%探测器与中心距离为r 时，这条探测器上从-1到1共2048个点， %转变为旋转角phi 下的坐标%将积分区间上的点转化为二维图像中的点%把坐标轴上的点的坐标转化为图像上的位置ϕϕϕϕϕϕϕϕcos sin )sin (cos )cos sin ()sin (cos t r y t r x t r l t r l +=-+=+-++=Φ+Φ=⊥tx = 256*(sqrt(2)*x + 1) + 0.5; ty = 256*(1 - sqrt(2)*y) + 0.5;5122*5.0225122*)1(5122*5.0225122*)1(-+--=+--=dy y dx x fx = floor(tx);fy = floor(ty);%把矩阵上的信息映射到图像上 ϕ%主要用差值法：%非边界点时用差值法，对于一个像素里的点用4的顶点的值去近似此点的值%if (tx<1 || ty<1 || tx>width || ty>height)mu = 0;elseif (tx == width && ty ~= height)mu = (fy+1 - ty)*(fx+1 - tx)*A(fx,fy) + (ty - fy)*(fx+1 - tx)*A(fx,fy+1);elseif (tx ~= width && ty == height)mu = (fy+1 - ty)*((fx+1 - tx)*A(fx,fy) + (tx - fx)*A(fx+1,fy));elseif (tx == width && ty == height)mu = A(fx,fy);elsemu1 = (fx+1 - tx)*A(fx,fy) + (tx - fx)*A(fx+1,fy);mu2 = (fx+1 - tx)*A(fx,fy+1) + (tx - fx)*A(fx+1,fy+1);mu = (fy+1 - ty)*mu1 + (ty - fy)*mu2;endint = int + dt*mu;End%设定足够长的滤波函数%P(i,:)*filter的长度分别为bins和2*bins-1；卷积的结果为的长度为3*bins-2,。

选取其中% [bins,2*bins]的范围内的值子函数2：function integ=RL_filter(P,view,bins,cellsize)filter_length = 2 * bins - 1;for i = 0:filter_length-1tmp = i - (filter_length - 1) / 2;if mod(tmp,2)==0filter(i+1) = 0.0;elsefilter(i+1) = -1.0 / (tmp * tmp * pi * pi) / cellsize;filter((filter_length - 1) / 2 + 1) = 0.25 / cellsize;endendfor i = 1 : viewP_filtered(i,:)=conv(P(i,:), filter);endinteg=P_filtered(:,bins:2*bins);子函数3：function Image=back_projection_circle(integ,cellsize,view,theta0,U)for height=1:512for width=1:512ee=0;for i=1:viewrr=(-U+(width-1)*cellsize+0.5*cellsize)*cos(i*theta0)+(U-(height-1)*cellsize-0.5*cellsize)*sin(i*theta0);%将图像上的点投影到接受器上，在接受器上的坐标（地址：有正有负）。