1、 平行束滤波反投影1100500121 赵伟伦准备知识:一维 Fourier 变换： dtetffFi2)()一维逆 Fourier 变换： efxf ti21)()()(且有： )(,11fFfFf 重要的性质：（卷积特性）；)()*(gff二维 Fourier 变换：；dXexffFxiR),(,21212 21),(),() 逆二维 Fourier 变换：；fxf xiR),(,21121 212),()(),( 中心切片定理：,)()(2fFfr其中 是 的 Radon 变换：,21x解释：一个二元函数的 Radon 变换关于 r 的一维 Fourier 变换与这个二元函数的二维Fou
2、rier 变换形式相等。滤波反投影：思路： )(),(121fFxfdrfFrdfef xrdefFfdeffFXxrrrirxirxirxixiR )(H,( fouriefourie),()(),( ,),(),(),(,),()(1),(20 2121 ),(,20 ),(,220),(,220 21),(,22121212121 变 化变 化 等 于 函 数 点 乘 后 的个 函 数 的 卷 积 的并 根 据 卷 积 的 性 质 ： 两设旋 转 角 为 为 坐 标 映 射 到 探 测 器 上， 设为用 极 坐 标 方 式 表 示 出 来（ 把 ， 可 知 ）， （ 由 于 中 心 切
3、片 定 理,),frG是滤波器(总结： drHf defFXfXXrir)(*,(,)0 20解释为：投影数据 ),(rf先进行滤波 )(*,(rHf在对滤波数据进行投影drfXr),0简单例子：（大圆与小圆）通过已得到的正投影round.dat经过滤波后，反投影后的图像。正投影数据：滤波图像：反投影后的图像：总的滤波反投影过程：1，得到图像的正投影2，滤波（投影与滤波器卷积）3，反投影（一个像素所有角度的滤波函数值的和再乘以 pi/views）viewsrHPdxfviews )(*,(,),021例子：GUI 界面实现：修改前：程序讲述：主程序：读入一张 512*512 的黑白图片，展示图
4、片；求正投，展示正投结果；求滤波后的结果并展示；求反投影后的的结果并展示。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
5、 = 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:文件ctzwl,s_png2) ;fp = fopen(s_png2,wb);fw
6、rite(fp,P,float64);fclose(fp);s_png = sprintf(angle.dat);%angles_png = strcat(D:文件ctzwl,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, imsh
7、ow(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
8、 = 0;for t = t1:dt:t2x = r*cos(phi) - t*sin(phi);y = r*sin(phi) + t*cos(phi);%探测器与中心距离为 r 时，这条探测器上从-1 到 1 共 2048 个点，%转变为旋转角 phi 下的坐标%将积分区间上的点转化为二维图像中的点%把坐标轴上的点的坐标转化为图像上的位置 cossin)i()cosin(tryxtltrtx = 256*(sqrt(2)*x + 1) + 0.5;ty = 256*(1 - sqrt(2)*y) + 0.5; 512*.0512*)(.dyxfx = floor(tx);fy = floor
9、(ty);%把矩阵上的信息映射到图像上%主要用差值法：%非边界点时用差值法，对于一个像素里的点用 4 的顶点的值去近似此点的值%if (txwidth | tyheight)mu = 0;elseif (tx = width elseif (tx = width elseif (tx = width 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;endi
10、nt = 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
11、) = -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;
12、for i=1:viewrr=(-U+(width-1)*cellsize+0.5*cellsize)*cos(i*theta0)+(U-(height-1)*cellsize-0.5*cellsize)*sin(i*theta0);%将图像上的点投影到接受器上，在接受器上的坐标（地址：有正有负） 。j=floor(rr/cellsize)+257;if(j511)ee=ee+0;elseaa=-U+(j-1)*cellsize+0.5*cellsize;%差值公式： ()()*cafcafbee=ee+(integ(i,j)+(integ(i,j+1)-integ(i,j)*(1/cellsize)*(rr-aa)*(pi/180);endendImage(width,height)=ee;endendend改一：改 t1,t2 和积分区域映射到二维图像的关系：让 t1=-sqrt（2） ，t2=sqrt（2）将正投影函数里： 把（-1 ，1）映射到第一个图像上。0.5; +y)-(1*256 =tx结果不对。改变反投影函数：结果
