用图像的投影数据进行重建程序

用图像的投影数据进行重建程序

宋利国赵向阳王春蔡国宝

摘要

本文通过引入Radon 变换来应用于CT图像重建问题，并且给出 Radon 变换在图像重建中的具体形式，并对滤波（卷积）逆投影法作了详细的研究，介绍其重建的基本思想和算法原理将问题转化为求解函数积分的形式。最后本文设计了一个人体脑部图像重建例子，通过matlab仿真说明如何投影才能重建准确的图像。

关键词：CT；图像重建；Radon变换；滤波逆投影法；matlab

1.问题重述

计算机断层成像技术（CT）是一种非介入式的检测技术，它极大地增强了人类观察物体内部结构的能力，在许多科学领域都得到了应用。特别在医学研究诊断中，它被用来作为一种获取人体内部信息的有效手段。我们在查阅许多资料，了解了CT成像的原理的基础上，选择采用在医学CT领域中的应用较为广泛，也是最基本最常用的图像重建算法──滤波逆投影法进行模型的仿真。

CT的工作原理就是投影重建（投影图像重建）。投影重建一般指从一个物体的多个（轴向）投影重建目标图像的过程。CT成像的基本数学原理是Radon 变换及其逆变换。目前，Radon变换及其逆变换是图像处理中的一种重要研究方法，许多图像重建便是有效地利用了这种方法，它不必知道图像内部的具体细节，仅利用图像的摄像值即可很好地反演出原图像。滤波逆投影法是当前用得较多的一种图像重建方法，在当代X 射线CT系统中几乎都用这种方法构成系统。它的特点是精度高，能快速实现。对于大量精确的投影数据来说，这是一种具有高效率的重建算法。滤波逆投影法又叫卷积逆投影法。这是因为频域上的滤波相当于空间域上的卷积运算。

我们通过引入Radon 变换来应用于CT图像重建问题，并且给出 Radon 变换在图像重建中的具体形式，对截面函数沿着特定直线进行积分就是它的 Radon 变换。滤波—逆投影法图像重建就是将截面函数沿若干个不同的角度下的特定直线进行积分产生的投影函数进行逆变换就得到了截面函数。滤波—反投影法能正确重建物体内部的吸收值图像，它把投影值按投影路线反过去赋予该路线上所有像元，使吸收值增加了该射线所经过的像元数目的倍数，经各个角度的投影反投回去与之叠加，最后能重建断面的图像。但由于反投影把投影路径的各处皆赋予该投影值，导致边缘较为模糊，所以通常把投影数值与某种校正函数卷积后再反投影，就能获得边缘清晰的图像。因为其中涉及到滤波函数的选取，也称为滤波反投影法。该重建方法兼顾了重建时间和重建质量两个方面，是医学上应用的最广泛的一种图像重建算法。

CT是X线照相术与复杂的计算机信号处理方法结合的产物，无论在医学放射诊断方面，还是在工业领域中均有着广泛的应用。采用滤波逆投影法成像技术，主要是因为医用CT可以采集到大量密集的投影数据，利用滤波逆投影法成像技术可以快速地得到具有一定质量的重建图像。

2.问题分析

2.1内容的选取

滤波逆投影法图像重建技术在医用CT 应用中的基本原理是由测量到的穿过人体横截面沿着许多直线的X 射线减的数据,重建出人体横截面的图像，是一种获取人体内部信息的有效手段，极大地增强了人类观察物体内部结构的能力，在医学成像方面发挥了巨大的作用。现代CT 成像的数学原理是Radon 变换及其逆变换。该变换是由函数在直线的线积分值来确定的,其逆变换就是由函数在空间所有直线上的线积分值确定这一函数（此函数对应实际中被扫描物体的密度函数或物体对X-射线的衰函数）。相对于早期的联立方程法和投影法克服了庞大的计算量和重建图像模糊精确度低的缺点。

2.2 影响因素的选取

在许多领域中，由于受客观条件的限制，经常会遇到不完全数据重建问题。

对于这类问题重建的讨论，不仅在理论上而且在实际应用中都有很重要的意义。这类问题可归结为以下几种情况：

[1]角度受限问题。这类不完全投影数据是由于在某些角度上无法采集到投影数据，从而导致投影角度覆盖范围缺损。在工业无损检测中，采集时间的限制或探测区域周围有障碍物都会导致这类不完全投影数据。

[2]外部问题。这类不完全投影数据是由于无法采集到穿过被探测区域某一部分的所有投影数据造成的。这种情况多发生在被探测物体内有对探测射线“不透明”物质存在时的情形。

[3]内部问题。这类不完全投影数据是由于射线的覆盖范围不能包含整个被探测的物体产生的。当被探测物体过大时或只想探测内部某一局部感兴趣区域时，就会发生这类不完全投影数据。内部问题也称为局部图像重建问题。在实际应用中，往往还会遇到以上三种类型中的组合情况，这时的重建问题更为复杂。对于不完全投影数据重建问题，基于各种重建判据的优化迭代重建算法是目前比较有效的种已知的先验知识进行外推，从而得到对缺失投影数据的估计，使迭代结果逐步向原图像逼近。

2.3 模型的选取

我们采用的是滤波反投影算法，该算法是目前广泛应用于所有直线透射断层成像的算法。此算法采用逆Radon变换以及投影定理，不仅极其准确而且保证了快速执行的稳定性。它是通过转换极坐标中Radon逆变换和将所有的定点依次积分得到了重建的图像来实现的。

3. 符号说明与基本假设

3.1 符号说明

由于变量与记号较多，为方便阅读，所有符号说明都将在其第一次出现时进行说明。

3.2 模型基本假设

H1：由于无法真正的利用医学仪器进行CT成像，我们采用的数据均是利用互联网搜集的，所以我们假设数据都是真实可靠的；

H2: 假设采集数据无不完全投影数据，可用于Radon变换及其逆变换；

H3：在无明显误差的情况下，我们假设重建后的图像与原图像基本吻合，建立的模型是正确的。

4. 模型处理的流程

反投影一般步骤为：

程序流程：

5.模型的建立与求解

5.1 Radon 变换及逆变换

设直角坐标系(,)x y 转动θ 角后得到旋转坐标系(,)x y ∧

∧

，由此得知

= x co s + ysin x θθ

∧

p (,)x θ∧

为原函数f (,)x y ∧

∧

的投影（f (,)x y 沿着旋转坐标系中x ∧

轴θ 方向的线积分）。 根据定义公式知其表达式为：

(,)(,)(co s sin ),,0d ef

p x f x y x y x d xd y x θδθθθπ∞

∞

∧

-∞-∞

=

+--∞∞≤≤??

这就是函数f (,)x y 的radon 变换。

从理论上讲图像重建过程就是逆radon 变换过程。即由投影函数p (,)x θ∧

逆向算出原函数f (,)x y 。其表达式为：

2

0(,)

1(,)(

)

2(cos sin )p x x

f x y d x d x y x

π

θθ

π

θθ∧

∧

∞

∧

∧

-∞

??=+-??

Radon 公式就是通过图像的大量线性积分来还原图像。为了达到准确的目的我们需要不同的θ建立很多旋转坐标系，从而可以得到大量的投影函数。为重建图像的精确度提供基础。 5.2滤波逆投影算法

有 radon 变换及投影定理可以方便的写出滤波逆投影方程：

(,)[co s(),][co s sin ](,)[co s()]f x y g r d g x y d p x h r x d x d πππθ?θθθθθθθ?θ

∧∧∧

+∞-∞

=

-=

+=

--?

?

??

方程中g 函数为角度为θ 时的累加函数。h 函数为滤波因子。p (,)x θ∧

为仪器得出的测量值函数。方程式中g 函数之起到一个中间变换的作用最终通过关系式

?

∞

∞

-=-)

,()(),('

''θθs g ds s s h s p

而被消掉。

由此我们就可以根据给定点(,)r θ在范围为(0,)π上进行积分。就可以得出该点位的函数值。将所有的定点依次积分就得到了重建的图像。 5.3滤波逆投影算法在matlab 中检测

在图像处理的工具箱中，MATLAB 提供了一个计算图像沿着指定方向上的投影的函数—radon 函数。iradon 函数可以实现radon 逆变换，radon 逆变换通常应用于X 线断层摄影术中，可以从投影数据中重构图像。

下面利用radon 函数和iradon 函数计算图像的投影并从投影中重建图像，Shepp-Logan 的大脑图作为测试图。函数radon 和函数iradon 的调用格式： [R,xp]=radon(I,theta) 计算图像I 在theta 向量所指定的方向上的radon 变换，I 表示待处理的图像，theta 表示radon 变换的方向角度，可以是标量或向量值，返回值R 的每一列对应图像I 在theta 某一角度的radon 变换值，xp 向量表示沿着x'轴对应的坐标值。IR =iradon(R,theta) 利用R 各列中投影值来构造图像I 的近似值。投影数越多，获得的图像越接近原始图像，角度theta 必须是固定增量的均匀向量。

6.仿真结果

6.1

仿真结果一

原始图

像

重建后图像

原图像 重建后图像

010*******-0.5

0.5投影函数(已补零)

010*******-1000

1000

S-L 卷积函数

020*******-100

100卷积结果

010*******

频谱

010*******

020*******

滤波逆投影法结果

6.2仿真结果二

原始图像

角度增值为10 时的iradon 变换图

像

角度增值为5时的iradon 变换图像

角度增值为2时的iradon 变换图像

θ

x '

经radon 变换后的图像

20

40

60

80

100

120

140

160

-150-100

-50

50

100

150

10

20

30

40

50

60

7.模型的分析评价

7.1模型的优点

引入了图像处理中的一种重要研究方法Radon 变换及其逆变换，这种方法不必知道图像内部的具体细节，仅利用图像的摄像值即可很好地反演出原图像；为了达到快速实现，使模型的精度高，我们采用了当前用得较多的一种图像重建方法滤波逆投影法。

7.2模型的不足与改进

由于时间、精力有限，我们的模型做得还不够精准，不够精湛，进一步的改进需要后续工作的努力。

8.总结

宋利国：

这次课程设计终于顺利完成了，在设计中遇到了很多编程问题，最后在老师的指导下，终于游逆而解。我们学得到很多实用的知识，在次我表示感谢！同时，对给过我帮助的所有同学和各位指导老师再次表示忠心的感谢！

课程设计是培养我们综合运用所学知识,发现,提出,分析和解决实际问题,锻炼实践能力的重要环节,是对我们实际工作能力的具体训练和考察过程.让我们进一步的了解了用图像的投影数据进行重建程序以及许许多多的知识内容。赵向阳：

本次课程设计颇费功夫，经过我们小组的通力合作，总算完成，从接到题目，分析课题，明确分工，查找资料，编写程序，合作的很默契，通过实验，图像处理的许多还没有学扎实的知识，通过本次试验，也很好的再温习了一遍，总之这个过程不是很顺利，学要大家的耐心才完成。

王春：

对于课程设计的顺利完成干到非常高兴，虽然在设计过程中遇到很多困难以及不懂的东西，但是通过小组成员的帮助，问题慢慢解答了，也重新认识到这门课程的各个方面的知识。在熟悉了书本知识的同时，也更加熟练了掌握了MATLAB相关技能。做完课程设计之后也了解到自己对于该门课程不足的地方，希望在以后的学习中更加努力，做到最好。

蔡国宝：

通过这次课程设计我对这门课程有了更多的了解，也知道自己不足的地方。在设计课程上还不熟悉，不够熟练说明自己的基础不是很好。通过这次课程设计学到了自己不足的地方，学会怎么样去克服遇到的困难。将课本上的理论知识变成实际操作，使自己的知识面更广，能力进一步提升。和同学的合作加强了自己的团队意识。总之，这次课程设计使我更加充分了解自己，使自己的各方面能力都得到提升。

参考文献：

[1]刘丹.计算机图像处理的数学和算法基础.国防工业出版社,2005.

[2]高欣.新型迭代图像重建算法的理论研究与实现.浙江大学,2004.

[3]张平.MATLAB基础与应用简明教程.北京航空航天大学出版社,2005.

[4]张含灵.MATLAB在图像处理中的应用.清华大学出版社,2008.

附录：

源程序代码一：

H=[0 0 0.92 0.69 90*pi/180 1;...

0 -0.0184 0.874 0.6624 90*pi/180 -0.98...

0.22 0 0.31 0.11 72*pi/180 -0.2...

-0.22 0 0.41 0.16 108*pi/180 -0.2;...

0 0.35 0.25 0.21 90*pi/180 0.1;...

0 0.1 0.046 0.046 0 0.2...

0 -0.1 0.046 0.046 0 0.2;...

-0.08 -0.605 0.046 0.023 0 0.1;...

0 -0.605 0.023 0.023 0 0.1;...

0.06 -0.605 0.046 0.023 90*pi/180 0.1];

angle=1;

N=100;

vstep=2/N;

ax=N/2+1;

for k=1:(180/angle+1)

theta=(k-1)*angle*pi/180;

for i=1:10

x0=H(i,1);y0=H(i,2);A=H(i,3);B=H(i,4);alpha=H(i,5);rho=H(i,6);

R=0;

forw=ax;

back=ax;

MM=sqrt(A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2-(R-x0*cos(theta)-y0*sin(t heta))^2);

NN=A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2;

g(i,ax)=rho*2*A*B*MM/NN;

for j=1:N/2

R=R+vstep;

forw=forw+1;

MM=sqrt(A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2-(R-x0*cos(theta)-y0*sin(t heta))^2);

NN=A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2;

g(i,forw)=rho*2*A*B*MM/NN;

R=-R;

back=back-1;

MM=sqrt(A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2-(R-x0*cos(theta)-y0*sin(t heta))^2);

NN=A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2;

g(i,back)=rho*2*A*B*MM/NN;

R=-R;

end

end

radon0(k,:)=real(sum(g));

end

%%%%%%%%%%%%%%%generateR-L function%%%%%%%%%%%%%%%% radon1=[zeros(k,N),radon0];

ax=N+1;

RL(ax)=1/(4*vstep^2);

forw=ax+1;

back=ax-1;

for k=1:N/2

n=2*k-1;

RL(forw)=-1/(n*pi*vstep)^2;

RL(back)=RL(forw);

RL(forw+1)=0;

RL(back-1)=0;

forw=forw+2;

back=back-2;

end

for k=1:(180/angle+1)

radon2(k,:)=conv(radon1(k,:),RL);

end

radonf=radon2(:,2*N:3*N);

%%%%%%%%%%%%%%%genrate S-L function%%%%%%%%%%%%

radon1=[zeros(k,N),radon0];

for v=1:(2*N+1)

n=v-N-1;

SL(v)=-2/(pi^2*vstep^2*(4*n^2-1));

end

for k=1:(180/angle+1)

radon2(k,:)=conv(radon1(k,:),SL);

end

radonf=radon2(:,2*N:3*N);

figure(1)

subplot(321)

plot(1:(2*N+1),radon1(1,:))

title('投影函数(已补零)')

subplot(323)

plot(1:(2*N+1),SL)

title('S-L卷积函数')

subplot(325)

plot(1:(4*N+1),radon2(1,:))

title('卷积结果')

Xradon1=fft(radon1(1,:));

subplot(322)

plot(1:(2*N+1),abs(Xradon1))

title('频谱')

XRL=fft(SL);

subplot(324)

plot(1:(2*N+1),abs(XRL))

Xradon2=fft(radon2(1,:));

subplot(326)

plot(1:(4*N+1),abs(Xradon2))

%%%%%%%%%%%%%%%%iradon%%%%%%%%%%%%%%%%%%%%

for k=1:(180/angle+1)

theta=(k-1)*angle*pi/180;

C=N/2-(N-1)*(cos(theta)+sin(theta))/2;

for i=1:N

R=(i-1)*cos(theta)+C;

n0=floor(R);

if n0>0&&n0<(N+1)

dot=R-n0;

I(i,l,k)=(l-dot)*radonf(k,n0)+dot*radonf(k,n+1);

else

I(i,l,k)=nan;

end

for j=2:N

R=R+sin(theta);

n0=floor(R);

if n0>0&&n0<(N+1)

dot=R-n0;

I(i,j,k)=(l-dot)*radonf(k,n0)+dot*radonf(k,n0+1);

else

I(i,j,k)=nan;

end

end

end

end

Ifinal=sum(I,3);

Gfinal=mat2gray(Ifinal);

Gfinal=imrotate(Gfinal,90);

figure(2)

imshow(Gfinal)

源程序代码二：

P=phantom(256); %用phantom 函数产生Sheep-Logan 的大脑图，n 为图像p 中的行列数，默认为256

imshow(P)

title('原始图像')

%以下为三种不同角度的投影模式

theta1=0:10:170;[R1,xp]=radon(P,theta1); %存在18个角度投影

theta2=0:5:175;[R2,xp]=radon(P,theta2); %存在36个角度投影

theta3=0:2:178;[R3,xp]=radon(P,theta3); %存在90个角度投影figure,imagesc(theta3,xp,R3);colormap(hot);colorbar;

%显示图像Sheep-Logan 的radon 变换

title('经radon变换后的图像')

xlabel('\theta');ylabel('x\prime'); %定义坐标轴

%用三种情况的逆radon 变换来重建图像

I1=iradon(R1,10);

I2=iradon(R2,5);

I3=iradon(R3,2);

figure,imshow(I1)

title('角度增值为10 时的iradon 变换图像')

figure,imshow(I2)

title('角度增值为5时的iradon 变换图像')

figure,imshow(I3)

title('角度增值为2时的iradon 变换图像')