扇束CT几何伪影的校正方法

《扇束CT几何伪影的校正方法》由会员分享，可在线阅读，更多相关《扇束CT几何伪影的校正方法（7页珍藏版）》请在装配图网上搜索。

1、第 5 期2011 年 5 月电子学报Vol . 39 No . 5May 2011ACTA EL ECTRONICA SINICA扇束 CT 几何伪影的校正方法王亮 ,张朋( 首都师范大学数学科学学院 ,北京 100048)摘要 : 实际工业 CT 系统很难满足经典 CT 图像重建公式所要求的射线源焦点 、转台中心及线阵探测器之间严格的对准关系. 不满足对准关系将导致经典 CT 图像重建公式所重建的图像中出现几何伪影. 本文引入可以描述一般情况下 CT 扫描几何关系的一组几何参数 ,给出了由金属丝的投影数据间接估计这些几何参数的方法 ,还给出了基 于这些几何参数的图像重建公式. 本文只假设转

2、台的旋转轴垂直于由射线源焦点与线阵探测器所形成的扇束 ,而不要 求扫描系统满足其它几何对准关系. 对于扫描系统不对准的情况 ,先按照本文方法估计几何参数 ,再将其代入本文提 出的图像重建公式 ,便可以正确重建出被测物的断层 CT 图像.关键词 :扇束 CT ; 图像重建算法 ; 扫描几何参数 ; 几何伪影校正中图分类号 : TP391文献标识码 :文章编号 :037222112 (2011) 0521143207ACo rre ctio n of Geo metric Artif a ct s fo r Fa n2B e a m CTWANG Liang ,ZHANG Peng( School

3、 of Mathematical Sciences , Capital Normal University , Beijing 100048 , China)Ab stract : For a real industrial CT system , it is very difficult to meet the rigorous aligned relations among the X2ray f ocus , the rotation center of turntable and the linear detector , which is required by the conven

4、tional image reconstruction f ormula . Misalign2 ment will result in geometric artif acts in the image reconstructed by the conventional image reconstruction f ormula . In this paper , we first introduce a set of geometric parameters which f ully describe general CT scanning configuration , and then

5、 p rop ose an indirect method to estimate the set of parameters f rom the p rojection data of a metal silk . We also p rop ose a generalized filtered back2p rojec2 tion ( FB P) reconstruction f ormula based on such set of parameters . In this paper , we assume that the rotation axis of the CT turnta

6、ble is vertical to the f an2beam f ormed by the X2ray f ocus and the linear detector array , but do not require that the CT scanning system satisfies other alignment relations . For a misaligned scanning system , we can first estimate the parameters by our method , and then reconstruct the image cor

7、rectly using our reconstruction f ormula with the estimated parameters .Key word s : f an2beam CT ; image reconstruction algorithm ; scanning geometric parameters ; geometric artif act correction参数 :射线源焦点到转台中心的距离、射线源焦点到探测器的距离. 此时 ,用理想的 FBP 重建公式所重建的 CT图像将包含几何伪影 ,导致图像分辨率降低 , 影响对图 像的正确判读. 当扫描视野与射线源焦点到探

8、测器的距 离之比较大时 ,图像的几何伪影将加重 , 特别是对图像中远离中心的部分. 因此 ,如何减轻或消除图像的几何伪影是 CT 设备研发的关键技术之一.有关减轻或消除 CT 图像几何伪影的研究结果公 开报道的相对较少. 1986 年 Gullberg 等2 提出一种带有转台中心偏移量、射线源焦点到转台中心和到探测器的距离等多个参数的扇束重建公式 ,并给出了参数的估计 方法. 该估计方法需要求解多个变量的非线性最优化问 题 ,鲁棒性较差. 陈明3 将转台中心的偏移归结为探测1 引言基于线阵探测器的断层工业 CT 目前仍然是工业 无损检测的主流设备. 三代扫描模式是断层 CT 最为常 用的扫描模

9、式. 三代 CT 扫描模式常用的图像重建公式 为经典的扇束滤波反投影( FBP) 重建公式. 该公式要求 CT 系统满足严格的几何成像对准关系1 : (1) 转台的旋 转轴垂直于由射线源焦点与线阵探测器所形成的扇束; (2) 射线源焦点与探测器中心的连线垂直于探测器 , 且 转台中心位于该垂线上. 此外 , 还要求精确测量 : ( 1) 射 线源焦点到转台中心的距离 ; ( 2) 射线源焦点到探测器 的距离. 然而实际 CT 系统机械安装时很难同时满足上 述对准关系 ,也难以直接精确测量图像重建所需的两个收稿日期 :2010211212 ;修回日期 :2011202218基金项目 :国家自然科

10、学基金(No . 60971131) ;北京市高层次人才项目(No . RHR20110509) ;北京自然科学基金(No . 3112006)器的偏转 ,给出了带有探测器偏转参量的重建公式 , 但没有考虑探测器中心的偏移量. 孙灵霞等4 设计了一 种射束中心位置测试装置 , 未涉及其它参数的测量. 对于锥束 CT 圆轨迹扫描系统 ,Noo 等人5 和陈炼等人6分别采用求质心的方法 , 求解扫描系统的部分几何参 数. 两篇文献都做了探测器平面平行于转轴的假定. 孙轴 , 其方向与 OS 相同; 过 O 且垂直于轴的直线为轴 , 使 O为右手直角坐标系. 转台坐标系 : 以 O 为原 点 , O

11、轴顺时针旋转一个角度得到 Ox 轴 , 建立右手直角坐标系 Oxy . 探测器坐标系 : 以 O 为原点 , 以探测器D所在直线为坐标轴建立坐标系 OD u .记 R = | OS | , O到 Ox 的顺时针夹角为, O到怡等人7 ,8O u 的顺时针夹角为. 设 O为 O与 O u 的交点 , 记在“4 小球模体所在平面垂直于转台平面”的DD假设下 ,利用一个角度下的投影数据给出了锥束 CT 扫描系统参数的一种测量方法. 针对锥束 CT 螺旋扫描系 统 ,牛小明等9 给出了带有探测器中心偏移量的重建方法. 该方法不需要对投影数据进行重排就能够减小 投影数据的偏移对重建图像的影响. 针对微焦

12、点锥束CT ,Patel 等10 提出一种不用测量模体的自校准方法 ,该方法对部分参数进行了测量与校正. 杨民等11 提出了 一种精确测量射线源焦点投影坐标的方法.针对三代 CT 图像的几何伪影问题 , 本文引入了可 以完全描述一般情况下的 CT 扫描几何关系的一组几何参数 ,给出了由金属丝的投影数据间接估计这组几何参数的方法 , 还给出了基于这组几何参数的图像重 建公式. 我们给出的几何参数间接估计方法的特点是:(1) 全部扫描几何参数均可以从金属丝模体的投影数 据中间接估计出来 , 克服了几何参数无法直接精确测量的困难 ; ( 2) 所有几何参数通过解析公式求出 , 不需 要求解非线性最优

13、化问题 , 因而求解结果具有很好的鲁棒性. 我们给出的基于这些几何参数的图像重建公式 ,只要求转台的旋转轴垂直于由射线源焦点与线阵 探测器所形成的扇束( 否则无法得到同一断层的扫描数据) ,而不要求扫描系统满足其它几何对准关系. 仿 真和实测数据的验证结果表明 : 在射线源焦点、转台中心和线阵探测器之间的几何位置不满足传统对准关系的情况下 ,利用我们的方法可以从金属丝的投影数据 中精确估计出几何参数 , 将这些几何参数代入我们给 出的图像重建公式就可以正确重建 CT 图像 ,有效消除 几何伪影. 当扫描视野与射线源焦点到探测器的距离 之比较大时 ,重建的 CT 图像同样具有较高的精度.2一般情

14、况下 CT 扫描几何参数的间接估计 方法本节首先给出一般情况下 , 如何选取几何参数和 坐标系来刻画 CT 扫描系统; 其次给出所选取的全部几 何参数的间接估计方法.211 一般情况下的扫描几何参数和坐标系如图 1 所示 , 以 S , O 和 OD 分别表示射线源焦点、 转台中心和探测器中心的位置. 我们将按如下方式建 立扇束坐标系 ( 固定坐标系) 、转台坐标系和探测器坐标系. 扇束坐标系 :选取 O 为原点; 过 O 与 S 的直线为h = | OOD| , D = | SO| ( 注意 : 当射线源焦点与转台中心的连线不垂直于线阵探测器时, D 表示转台中心在 探测器上投影点 O到射线

15、源焦点 S 的距离 , 而不是射 线源焦点到探测器的距离) . 固定参数 R , D , h ,与角 度参数一起 , 唯一刻画了 CT 扫描状态.212 任意点 P0 的投影地址设 P0 为物体上远离转台中心的任意一点, 其在Oxy 下的坐标记为( x0 , y0) ; 当 Ox 与 O夹角为时 , P0在 O下的坐标记为 (0 ,0 ,) ; P0 在 OD u 上的投影 地址记为 u0 ,. 下面用( x0 , y0) 及参数 R , D , h ,和来 表示 u0 ,和(0 ,0 ,) .由 S 发出且过点 P0 的射线交探测器于点 U0 , 在O坐标系下可表示为OU0 , = R (

16、0 ,1) + t0 ,( (0 ,0 ,) - R (0 , 1) )( 1)另一方面 , 点 U0 ,也可表示为OU0 , = - ( D - R) ( 0 ,1) + ( h + u0 ,) ( n1 , n2) ( 2) 其中( n1 , n2) = ( cos, sin) . 由式( 1) 与( 2) 求解可得 P0 在 OD u 上的投影地址 D0 ,u0 , = - h + n( 3)0 ,) 20 ,1 ( R -+ n0 , = x0cos+ y0 sin其中( 4)0 , = y0cos - x0 sin2138 个对称点的投影地址满足的方程组为了简单起见 , 我们做变量替

17、换. 令n1n2x0y0n1 = D , n2 = D , x0 = R , y0 =R( 5)0 , = x 0cos+ y0 sin0 , = y 0cos- x0 sin第 5 期王 亮 :扇束 CT 几何伪影的校正方法1145记( x0 , y0) 所在的坐标系为 Oxy , (0 ,0 ,) 所在的坐标系为 O. ( x 0 , y0 ) 与 ( x0 , y0) , (0 ,0 ,) 与 (0 ,0 ,) 的关系是固定比例的缩放关系 , 图像经过该缩放后结构不变. 原来的几何参数 R , D , h , n1 , n2 中 n1 与 n2 是不独立的 , 而新的几何参数 R , h

18、 , n1 , n2 是彼此独 立的. 由 R , h , n1 , n2 与角度参数一起 , 可以唯一刻画 CT 扫描状态. 做变量替换后 , 由式 ( 3) 可以得到 (0 , 0 ,) 在探测器坐标下的投影地址0 ,u3 + h2 - m1 + m2u7 + h = -(13)2 + m - m1 2这 4 个方程中包含 3 个未知数 h , m , m , 从中任选 3 个1 2方程联立 , 便可消去 m1 和 m2 求得 h . 于是 h 可有以下4 个表达式h =- ( u0 + u4) U1 U3 + 2 U0 ( u3 u5 -u1 u7)- 12 U U + 2 U U -

19、U0 (3)(14)1 31h =- ( u2 + u6) U1 U3 - 2 U2 ( u5 u7 -u1 u3)- 12 U1 U3 - 2 U2 ( U1 + U3)u0 , = - h +( 6)(15)n1 ( 1 - 0 ,) + n20 ,h =- ( u1 + u5) U0 U2 + 2 U1 ( u0 u2 -u4 u6)当 Ox 与 O 夹角为时 , ( x 0 , y0 ) 在 O下的坐标为(0 ,0 ,) . 令 ( x0 , y 0 ) 顺时针每次旋转 / 4 , 旋转 7- 12 U0 U2 - 2 U1 ( U0 + U2)(16)h =次后 得 到 8 个 对

20、称 点. 该 8 个 点 相 应 的 (0 ,+ k ,- ( u3 + u7) U0 U2 - 2 U3 ( u2 u4 -u0 u6)4- 12 U0 U2 - 2 U3 ( U0 - U2)0 ,+ k) 之间有如下关系(17)4其中 U0 = u0 - u4 , U1 = u1 - u5 , U2 = u2 - u6 , U3 = u3 -u7 . 由于投影地址之间存在着内在的相容关系, 因此在cos k- sin k44(0 ,+ k ,0 ,+ k) = (0 ,0 ,),sin kcos k44( x , y 的投影地址 u ( k = 0 ,1 ,2 ,7) 不含误差时 , 由

21、0)0k44, 7式( 1417) 计算出的 h 是相同的. 在投影地址 u 存在kk = 0 ,1 ,2 ,( 7)误差时 , 将式( 1417) 计算出的 h 进行平均 , 可以降低投影地址 uk 误差引起的 h 的计算误差.求解出参数 h 之后 , 现在考虑求解参数 n1 和 n2 . 为此 , 联立 u0 与 u4 , u2 与 u6 分别解出两组0 ,和0 , 的表达式利用式(7) , 上述 8 个对称点的投影地址可表示为0 ,+ k4 ,u = - h + 0 ,+ kn ( 1 - ) + n 0 ,+ k 2 0 ,+ k144k = 0 ,1 ,2 ,4, 7( 8)214由

22、投影地址估计几何参数的公式本小节将导出由 8 个对称点的投影地址估计几何 参数 h , n1 , n2 和 R 的公式.- 10 , = - 2 H0 H4 n1 U0(18)0 , =0 , =- 1- 2 H0 H4 n2 + H0 + H4U0- 1为简单起见 , 令 uk = u , k = 0 ,1 ,2 , 7 . 由 u00 ,+ k2 H2 H6 n2 - H2 - H6U24(19)与 u4 , u1 与 u5 , u2 与 u6 , u3 与 u7 可得如下方程组- 10 , = - 2 H2 H6 n1 U2其中 H0 = h + u0 , H1 = h + u1 , H

23、2 = h + u2 , H3 = h + u3 ,H4 = h + u4 , H5 = h + u5 , H6 = h + u6 , H7 = h + u7 . 从中消去0 ,和0 , 即可求解出 n1 和 n2u + h0 ,+ k 4 + hu 0 ,+ ( k + 4) 4n1 ( 1 - 0 ,+) + n20 ,+n1 = U0 U2 ( H H ( H + H ) - H H ( H + H )0 ,+ k( k + 4) ( k + 4) 0 4 2 6 26 0 4444(2 H02 H42 U22 + 2 H22 H62 U02) - 1=(20),n1 ( 1 - k)

24、+ n2 kn2 = ( U02 H H H + H+ U 2 H H H + H0 ,+ 40 ,+ 46 (6)4 (4) )0 ,+ ( k + 4) 2 22 0 04(2 H02 H42 U22 + 2 H22 H62 U02) - 1(21)( 9)k = 0 ,1 ,2 , 3由另外 4 个角度下的投影 u 与 u , u 与 u 同样可以1 5 37n2n2解出两组0 ,和0 ,的表达式 , 进而可解得 n1 和 n2令 m1 = n 0 , +0 ,m2 = n 0 , - 0 ,11n1 = U1 U3H1 H5 ( H3 + H7) - H3 H7 ( H1 + H5)则

25、方程组( 9) 可化简为2 H12 H52 U32 + 2 H32 H72 U12- 1(22) u0 + h1 - m2u4 + h = - 1 + m2( 10)U12 H3 H7 ( H3 + H7) + U32 H1 H5 ( H1 + H5)n =22 2 22 2 2 - 1 u2 + h1 - m12 H1 H5 U3 + 2 H3 H7 U1(23)u6 + h = - 1 + m1( 11)同样 , 由于投影地址间存在着内在相容关系 , 在投影地址 uk ( k = 0 ,1 , 2 ,7) 不含误差时 , 由式( 20) ( 22) , u1 + h 2 - m1 - m2

26、( 12)u5 + h = -( 21) ( 23) 分别计算出的 n , n 是相同的. 在投影地址 u2 + m + m1 21 2k存在误差时 , 将式 ( 20) ( 22) , ( 21) ( 23) 计算出的 n1 , n2进行平均 , 可以降低投影地址 uk 误差引起的 n1 , n2 的 计算误差.在实际应用中 , 我们选取一根细的铅或钨的金属 丝( 或水银温度计) 尽量垂直地立于转台上远离中心的位置 , 并保持金属丝 ( 或水银温度计) 的投影位于探测器视野内. 我们采用三代 CT“转停”采样模式 , 即转台旋 转一个角度后停止 , 探测器开始采集穿过被测物的射 线 , 当采

27、集结束后再旋转到下一个角度位置采集数据,直至完成一周采样. 旋转一周的角度采样数为8 的倍数 ( 如 720 或 1800) . 从金属丝( 或水银温度计) 扫描数据图 像( Sinogram) 中提取出其中心的投影地址. 利用式( 1423) , 从每个 0 ,/ 4) 相应的 8 个投影地址 , 分别计 算出参数 h , n1 和 n2 ; 为提高数值计算的稳定性 , 将相 应不同求得的参数 h , n1 和 n2 进行平均.接下来考虑计算参数 R . 注意到由 h , n1 , n2 和式(18 ,19) 可以求出0 ,和0 , 因此 , 利用两根标定了距 离的金属丝就可以求出参数 R

28、, 具体方法如下 :设两根平行的金属丝之间的距离为 d , ( x1 , y1 ) 和( x2 , y2) 分别为两根金属丝中心在 CT 断层中的坐标 , 它 们在 O下的坐标分别记为(1 ,1 ,) 和 (2 ,2 ,) ,利用 213 小节中变换式( 5) 可得 , 在 O下的坐标分别记为(1 ,1 ,) 和(2 ,2 ,) . 易见有下列关系下面通过建立平行束下参数( r , ) 与扇束下参数( u ,) 之间的对应关系 , 并利用坐标变换 , 推导基于参 数 h , n1 , n2 的扇束 CT 滤波反投影图像重建公式.第一步 , 推导 ( r , 0的方法可得. 用求极值R ( h

29、+ u) n1(28)r =( h + u) 2 n12 +2D - ( h + u) n2设 Q 为原点 O 到射线 SU 的垂足 , 可求得 OQ 在 O下的表达式为D - ( h + u) n2R ( h + u) n1(,) =,2( h + u) 2 n12 +D - ( h + u) n22D - ( h + u) n2RR -(29)2( h + u) 2 n12 +D - ( h + u) n2另一方面cossin- sin cosr ( cos , sin ) = (,)由式( 2830) 解得(30)( h + u) n1cos+ ( D - ( h + u) n2) si

30、n = arctan(31)- ( h + u) n1 sin+ ( D - ( h + u) n2) cos式( 28) ( 31) 即是参数( r , ) 与参数( u ,) 之间的关系.第二步 , 确定 d rd 与 d ud的关系. 由式( 28) 和式( 31) 可求得雅克比行列式( x1 , y1) - ( x2 , y2)d =(1 ,1 ,) - (2 ,2 ,)(1 ,1 ,) - (2 ,2 ,)= R( 24)D - ( h + u) n25 r , RDn1 ( )利用本节方法可以求出几何参数 h , n1 , n2 , 而由 h , n1 ,n2 和式( 18 、19

31、) 可以求出(1 ,1 ,) 和(2 ,2 ,) , 因此可得到射线源焦点到转台中心的距离为d= (32)5 ( u ,)322 n 2 +D - ( h + u) n2( h + u)21第三步 , 将 x0cos + y0 sin - r 用 u , R , D , h , n1 ,n2 表示. 点 P0 在 Oxy 下的坐标为 ( x0 , y0 ) , 由式 ( 4) 可得(0 ,0 ,) 所表示的( x0 , y0) 为x0 = 0 ,cos- 0 ,sinR =( 25)(1 ,1 ,) - (2 ,2 ,)(33)3 基于参数 h , n1 , n2 的图像重建公式下面我们从平行

32、束的滤波反投影 CT 图像重建公 式出发 , 推导一般情况下扇束扫描基于参数 h , n1 , n2 的滤波反投影图像重建公式.平行束的滤波反投影 CT 图像重建公式可以写成 如下形式y0 = 0 ,sin+ 0 ,cos由式( 2830) 可得( cos , sin ) = - ( h + u) n1 sin+ ( D - ( h + u) n2) cos,( h + u) 2 n12 +2D - ( h + u) n2( h + u) n1cos+ ( D - ( h + u) n2) sin(34)2 + ( h + u) 2 n12 + ( D - ( h + u) n2) 2 1f

33、( x0 , y0) = 2 f ( r , ) H( x0cos + y0 sin - r) d rd 0 - 由式( 28 , 33 ,34 得)( 26) 其中 , ( x0 , y0) 为待重建点 , r 为坐标原点 O 到射线 SU 的距离 , 为射线法向与 x 轴方向的夹角 ( 如图 1 所- hn1 ( R - 0 ,) + ( D -hn2) 0 ,x0cos + y0 sin - r =n1 ( R - 0 ,) + n20 , u- 12 22 ( h + u)n1 + ( D - ( h + u) n2)(35)2+ 示) , f ( r , ) 是以( r , ) 为变

34、量的投影数据 , H( r) = 由投影地址公式( 3) 得- hn1 ( R - 0 ,) + ( D - | | ei2rd.hn2) 0 ,OP第 5 期王亮 :扇束 CT 几何伪影的校正方法1147n1 ( R - 0 ,) + n20 , u0 ,所代表的长度进行标定时才需要 R 的值. 如果模体的长度已知 , 通过比较该模体 CT 图像的像素单位长度和 它的实际长度 , 也可推出 R 的值.4实验验证本节我们分别采用仿真数据和实测工业 CT 数据 ,验证本文提出的几何伪影校正方法的有效性.411仿真数据验证用直径为 0175mm 的小圆模拟金属丝的横截面. 系 统的 扫 描 几 何

35、 参 数 分 别 设 定 为 R = 1000mm , D =1200mm. 线阵探测器由 1400 个长度为 0125mm 的单元组成 ,总长度 350mm. 此时放大比为 112 , 像素单元尺寸为012mm ,扫描视野直径为 2881614mm. 角度采样数为 1800 .为验证本文参数间接估计方法的正确性和稳定 性 , 选取 3 组参数进行数值实验. No . 1 组参数值及金属丝的 位 置 分 别 设 定 为 : h = 2mm , n1 = 0199996 , n2 =0100873 (= 015) , x0 = 130mm , y0 = 40mm , 模拟金属丝 的断层扫描数据

36、, 如图 2 所示. 图中横向为探测器坐标 ,纵向为采样角度. No . 2 组和 No . 3 组参数值及金属丝的 位置分别设定为 : h = 4mm , n1 = 0199985 , n2 = 0101745 (= 1) , x0 = - 120mm , y0 = - 60mm 和 h = 6mm , n1 =( 36)将式(36) 代入(35) 右端 , 则n1 ( R - 0 ,) + n20 ,( u0 , - u)x0cos + y0 sin - r =2( h + u) 2 n12 +D - ( h + u) n2( 37)又滤波函数 H( 3 ) 为负二次齐次 , 因此有H(

37、x0cos + y0 sin - r) =2( h + u) 2 n12 +D - ( h + u) n2H( u0 , - u)( 38)n1 ( R - 0 ,) + n20 ,2将式( 32 , 38) 代入平行束重建公式( 26) , 得到带有扫描几何参数 R , D , h , n1 , n2 的扇束滤波反投影重建 公式2 + 1f ( x0 , y0) = 2 p ( u ,)0 - D - ( h + u) n22( h + u) 2 n12 +D -( h + u) n2RDn1 H( u0 , - u)d ud( 39)n1 ( R - 0 ,) + n20 ,2其中 p (

38、 u ,) 是以( u ,) 为变量的扇束投影数据.将重建公式( 39) 右端第一个分式的分子分母同时 除以 D , 第二个分式的分子分母同时除以 R2 D2 , 公式两 端再同时乘 R , 由变量替换式( 5) 可得0199939 , n = 0103489 (= 2) , x = 95mm , y= - 95mm.200其断层扫描数据图与图 2 相似 , 只是起始位置不同.2 1Rf ( x0 , y0) = 2 p ( u ,)0 - 1 - ( h + u) n2( h + u) 2 n12 + ( 1 - ( h + u) n2) 2 n1 H( u0 , - u) 2 d ud(

39、40)( n1 ( 1 - 0 ,) + n20 ,)将式( 40) 右端记作新的函数 f ( x0 , y 0 ) , 得到基于参数 h , n1 和 n2 的扇束 CT 滤波反投影图像重建公式2 f ( x0 , y 0) = 12 0 - p ( u ,)1 - ( h + u) n2采用图像处理的方法提取质心 , 可以得到投影图中金属丝中心的投影地址 , 利用本文第 2 节提出的间接( h + u) 2 n12 + ( 1 -n1 H( u0 , - u)( h + u) n2) 2估计方法计算出参数 h , n 和 n . 为保证计算结果的稳1 2定性 , 对于每组实验 , 将 36

40、0 度内不同角度对应的估计值进行平均后作为参数的最终计算结果, 结果见表 1 .由表 1 中 No . 1 组参数值的实验结果可知 : (1) 计算 出的 h 平均值为 211165mm , 与设定值 2mm 相比 , 误差为011165mm. 实验所设定的线阵探测器单元长度为0125mm , 该误差相当于探测器单元长度的 01466 倍. ( 2)根据估计值 n1 = 813328 10 - 4 , n2 = 718246 10 - 6 可以2 d ud( 41)( n1 ( 1 - 0 ,) + n20 ,)于是有如下关系f ( x0 , y0) = Rf ( x0 , y0)( 42)从

41、关系 x0 = x0/ R , y0 = y0/ R 和式 ( 42) 可知 , R 的值不影响f ( x 0 , y 0) 和 f ( x0 , y0) 量化后的图像值 , 而只影响它 们相应图像的整体缩放比例 , 即 f ( x0 , y0 ) 和 f ( x0 , y0 )相应图像是相似的 , 因此可将重建 f ( x0 , y0 ) 改为重建f ( x0 , y 0) .与式(39) 相比 , 式 ( 41) 中只含有 h , n1 和 n23 个参 数 , 而且这 3 个参数是彼此独立的. 此外 , 从式 ( 41) 可 知 , 在图像重建中并不需要知道 R 的值. 对图像上像素1计

42、算 出 D = ( n 2 + n 2 ) - = 12001024mm , 与 设 定 值21 21200mm 相比 , 误差为 01024mm. ( 3) 根据 n , n 可以计1 2算出探测器的偏转角 = arctan ( n2/ n1 ) = 010094 弧度(01538 度) , 与设定值 015 度相比 , 误差为 01038 度. 类 似地 , 可以计算出表1 中No . 2 ,No . 3组参数值所对应的h , 根据 n1 , n2 计算出 D 以及探测器的偏转角.表 1 参数间接估计结果通过本文第 2 节方法计算出的 h , n1 和 n2 , 如表 21所示. 由参数

43、n 和 n 可以计算出 D = ( n 2 + n 2) - =21 21 2No . 参数设定值估计值( 平均值)绝对误差10191165mm , 探测器偏转角= arctan ( n / n ) = 01024852 1h ( mm)n 1n 228 . 333 10 - 47 . 272 10 - 62 . 11658. 3328 10 - 47. 8246 10 - 60 . 11652 . 16719 10 - 85 . 52508 10 - 7弧度( 11424 度) .为验证参数估计结 果的正确性 , 在 同 样 的扫描条件下 , 扫 描 两 个1表 2 CT 系统参数估计结果h

44、 ( mm)n 1n 248 . 332 10 - 41 . 454 10 - 54. 127238 . 33178 10 - 41. 4356 10 - 50 . 127238 . 21602 10 - 81 . 88642 10 - 72易拉罐 一个靠近转台(h ( mm)n 1n 268 . 328 10 - 42 . 908 10 - 56. 128488 . 32762 10 - 42. 9267 10 - 50 . 12848- 6. 34234 10 - 81 . 84054 10 - 7中心 , 一个远离 转 台 中3心) , 获得其扫描数据. 图 3 给出了下列三种方法的重建

45、结果 :方法 1 为理想的 FBP 重建算法; 方法 2 为调节 探测器中心偏移量(offset) 的重建算法; 方法 3 为本文提412实测数据验证实际 CT 系统所用的探测器单元宽度为 01083mm ,探测器单元个数为 3710 ; X 射线源焦点的 1/ 10 高宽为3mm ;射线源焦点到转台中心的距离约为 800mm ; 射线 源焦 点 到 探 测 器 的 距 离 约 为 1000mm. 将 直 径 约 为015mm 的细铅丝立于转台上远离旋转中心的位置处进行扫描 ,角度采样数为 1800 . 从铅丝的扫描数据中可以 提取出不同角度下的中心位置.出的基于参数 h , n 和 n 的扇

46、束滤波反投影重建算1 2法. 图 3 ( a) 为方法 1 的重建结果 , 从其局部放大图 3( a1) 、3 ( a2) 中可以看出易拉罐壁存在明显的几何错位 伪影 , 这是由扫描几何误差所引起的. 图 3 ( b) 为方法 2的重建结果 , 从其局部放大图 3 ( b1) 、3 ( b2) 中可以看出靠近转台中心处的易拉罐壁基本正确, 但远离转台中 心处的易拉罐壁仍有明显的几何错位伪影. 图 3( c ) 为方法 3 的重建结 果 , 从 其 局 部 放 大 图 3( c1) 、3 ( c2) 中可以看出 重建的易拉罐壁在靠近转台中心处和远离转台 中心处均无明显伪影 , 说明本文方法可以对

47、扫描 系统几何位置偏差所引起的伪影进行有效的校正.5总结本文提出了扇束 CT 几何伪影的校正方法 ,该 方法只要求转台的旋转 轴垂直于由射线源焦点 与线阵探测器所形成的 扇束 ,而不要求扫描系统 满足其它几何对准关系. 本文重建公式所需 的几何参数可以从金属 丝的投影数据中间接估 计 ,从而克服了经典重建 公式中所需参数无法直 接精确测量的困难. 由于参数估计值( 平均值)h (mm)0. 328265n19 . 809 10 - 4n22 . 438 10 - 5第 5 期王 亮 :扇束 CT 几何伪影的校正方法1149所有参数可通过解析公式求出 , 不需要求解非线性最优化问题 ,因而求解结

48、果具有很好的鲁棒性. 而且该间 接估计方法还降低了对模体制作精度和复杂度的要求.实验表明在射线源焦点、转台中心和线阵探测器不 满足传统对准关系的情况下 ,利用本文提出的几何参数的间接估计方法及基于这些几何参数的 CT 图像重建公 式 ,仍可正确重建 CT 图像. 当扫描视野与射线源焦点到探测器的距离之比较大时 ,重建的 CT 图像也具有较高 精度.本文间接估计方法假定了角度采样数是8 的倍数 ,即投影具有 8 对称性. 实际系统如采用闭环控制的“转 停”模式 ,则可以精确做到角度采样数是8 的倍数. 如果采用“连续”模式 , 转台的运动和探测器采样可能不耦 合 ,造成一周内的角度采样数可能不是

49、8 的倍数. 此时 ,可以先采用闭环控制的“转停”采样模式采集金属丝的 扫描数据 ,进而间接估计出全部几何参数. 在“连续”扫描模式下 ,将这些几何参数代入重建公式中即可重建高精度的 CT 图像. 也可以将 ( 8) 改为具有准对称关系 的 8 个投影 ,导出一般的几何参数间接估计公式.致谢 :感谢首都师范大学检测成像实验室张慧滔、 陈德峰老师以及朱溢佞、赵云松同学的有益讨论及在 数值实验方面给予的帮助 , 感谢赵星老师对本文修改 的建议. 6 陈炼 , 吴志芳 , 刘锡明 , 等. 锥束 CT 系统几何参数校正的解析计算J . 清华大学学报 ( 自然科学版) , 2010 , 50 ( 3)

50、 :418 - 421.Chen Lian , Wu Zhif ang ,Liu Ximing , et al . Analytical geometric parameter calibration algorithm f or cone2beam CT J . J ournalof Tsinghua University ( Science and Technology ) , 2010 , 50(3) :418 - 421. (in Chinese) 7 Sun Y , Hou Y , Zhao F Y , et al . A calibration method f or mis2 a

51、ligned scanner geometry in cone2beam comp uted tomograp hy J . ND T & E International , 2006 , 39 (6) :499 - 513. 8 Sun Y , Hou Y , Hu J S . Reduction of artif acts induced by mis2 aligned geometry in cone2beam CT J . I EEE Transactions on Biomedical Engineering , 2007 , 54 (8) :1461 - 1471. 9 牛小明 ,

52、 曾理 , 刘宝东 , 等. 偏心螺旋锥束 CT 的 Katsevich修正重建算法J . 计算机工程 , 2009 , 35 (5) :215 - 217.Niu Xiaoming , Zeng Li , Liu Baodong , et al . Katsevich correc2 tion and rectification algorithm f or off2center spiral cone2beam CTJ . Comp uter Engineering , 2009 , 35 ( 5 ) : 215 - 217. ( in Chinese)10 Patel V , Chity

53、ala R N , Hoff mann K R , et al . Self2calibration of a cone2beam micro2CT system J . Med Phys , 2009 , 36 (1) :48 - 58. 11 杨民 , 刘静华 , 李保磊等. 三维显微 CT 扫描系统的 X 射线源焦点投影坐标测量方法J . 光学学报 , 2009 , 29 ( 5) :1275 - 1280.Yang Min ,Liu J inghua ,Li Baolei , et al . Measurement method of X2ray f ocus p rojection

54、coordinates of three2dimensional mi2 cro2comp uterized tomograp hy scanning system J . Acta Opti 2ca Sinica , 2009 , 29 (5) :1275 - 1280. (in Chinese)参考文献 1 张朝宗 , 郭志平 , 张朋 , 等. 工业 CT 技术和原理 M .京 :科学出版社 , 2009. 122 - 135.北作者简介 2 Gullberg G T , Crawf ord C R , Benjamin M W . Reconstructionalgorithm f o

55、r f an beam with a displaced center2of2rotation J . I EEE Trans Med Imag , 1986 , 5 (1) :23 - 29. 3 陈明. 工业 CT 视野拓展方法及图像伪影校正研究 D .北京 :首都师范大学 , 2008. 6Chen Ming , Research of field2of2view extension and image arti2 f acts correction in industrial comp uted tomograp hy D . Bei 2 jing : Capital Normal U

56、niversity , 2008. (in Chinese) 4 孙灵霞 , 李炬 , 周江 , 等. 工业 CT 系统中射线源的精确定位J . 核电子学与探测技术 , 2007 , 27 (4) :716 - 718.Sun Lingxia ,Li J u , Zhou J iang , et al . Precise p ositioning of X2 ray source in industrial CT system J . Nuclear Electronics & Detection Technology , 2007 , 27 (4) :716 - 718. (in Chinese) 5 Noo F , Clackdoyle R , Mennessier C. Analytic method based onidentification of ellipse parameters f or scanner calibration in cone2beam tomograp hy J . Phys Med Biol , 2000 , 45 : 3489 -3508.王 亮