滤波反投影

平行束滤波反投影

1100500121 赵伟伦 准备知识:

一维Fourier 变换：

dt e

t 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 --⋅=⋅=

重要的性质：（卷积特性）

)(~)(~)*(ωωg

f g f F ⋅=； )(~)(~)(ωωg

f g f F *=⋅ 二维Fourier 变换： dX e x x f f f F x x i R ),(),(22

121221212),(),(~)(⋅-⎰=

=ωωπωω； 逆二维Fourier 变换： Ω==⋅-⎰d e f f F x x f x x i R ),(),(22

1122121212),(~)~(),(ωωπωω； 中心切片定理：

),)(ˆ()(2ϕωωf

F f F r =Φ, 其中),(ˆϕr f 是),(2

1x x f 的Radon 变换： 解释：一个二元函数的Radon 变换关于r 的一维Fourier 变换与这个二元函数的二维Fourier 变换形式相等。

滤波反投影：

思路：

)(),(121f F F x x f ⋅=-

()

()[][]ϕ

ϕωωϕ

ωϕωϕωωϕωϕωϕωωωϕωωϕωϕωωϕωϕωωωϕωωωππωωππωωππωωππωωπd r f F r d f

F F d d e f

F x x r d d e f

F 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 r

x 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 f

d d

e

f 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

πϕϕ

ϕϕπ

∆=*=∑⎰=)(*),()(),(),(00

21

例子：

GUI界面实现：

修改前：

程序讲述：

主程序：

读入一张512*512的黑白图片，展示图片；求正投，展示正投结果；求滤波后的结果并展示；求反投影后的的结果并展示。

view=360;

bins=512;

U=1.5; %U=view_radius

P=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:view

phi = (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;

end

end

outtime = 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');%angle

s_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,[])