Advances in Applied Mathematics
Vol.06 No.04(2017), Article ID:21395,8 pages
10.12677/AAM.2017.64059

Formulas of the Length and Index for the Intersection of a Line and a Fixed Grid in Three-Dimensional Space

Xiaoyuan Zhang, Caifang Wang, Xudao Yin, Zhaoliang Xu

College of Arts and Sciences, Shanghai Maritime University, Shanghai

Received: Jun. 23rd, 2017; accepted: Jul. 15th, 2017; published: Jul. 18th, 2017

ABSTRACT

In this paper, we introduce the principle of imaging for X-ray CT and try to convert the reconstruction problem to a linear system Ax=b. We emphasize on the construction of A. Using a reconstruction example, we provide formulas of length and index of the intersection of a line and a fixed grid in three-dimensional space. Finally, we try to verify the formulas with MATLAB programming.

Keywords:Digital Image Reconstruction, X-Ray CT, Grid in 3D, Index

直线与三维固定网格的交线长度与 索引公式计算

张晓媛,王彩芳,殷绪导,徐兆亮

上海海事大学文理学院,上海

收稿日期:2017年6月23日;录用日期:2017年7月15日;发布日期:2017年7月18日

摘 要

本文首先介绍了数字图像重建中最典型的计算机断层成像技术(X-ray Computerized Tomography)的成像原理,并考虑将重建问题转化成线性方程组Ax = b的求解,重点讨论了系数矩阵A的生成。以三维实际重建问题为例,本文给出了计算直线与三维固定网格的交线长度与索引公式。最后利用MATLAB实现这一计算过程,并分析算法运行的效率。

关键词 :数字图像重建,X-Ray CT,三维固定网格,索引

Copyright © 2017 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

数字图像重建是指根据物体边界测量数据,获取物体内部结构的一种方法 [1] 。常见的成像技术有全息成像,电子显微,数码相机,核磁共振,超声波,计算机断层层析成像技术等等。图像重建最典型的应用是医学上的计算机断层层析成像技术 [2] ,这是一种电子计算机技术与X射线检查技术相结合的产物,因其扫描时间快,图像清晰等特点,在成像领域有着非常重要的作用。

X射线在物体内部沿直线传播,衰减满足Beer定理 [3]

其中为物体对X射线的衰减系数,为入射光强度,为出射光强度。X-Ray CT利用这一原理,获取多组入射光强和出射光强度,得到一系列的投影数据,并根据投影数据值重建物体内部的吸收系数.

常见的X-Ray CT重建算法有迭代算法 [4] ,分析算法(滤波反投影 [5] ,反投影滤波 [6] )等。随着计算机技术的高速发展,迭代重建算法因其利用计算机运算速度快,适合做重复性操作的特点被人们更广泛地应用在图像重建中。迭代算法首先涉及到问题的离散化,以三维图像为例,当待测物体所在的区域被离散成三维固定网格后(共个像素),如图1所示,假设物体在离散后每个格点上的吸收系数是一个常数,沿着第条射线上的投影值,可表示成,其中为第条射线与第个像素的交线长度,这样就将重建问题转化成的求解 [7] [8] 。我们的任务是确定矩阵的第个行向量中的非零元以及非零元所在列的索引。

作者吴大瑞,何钦铭在 [9] 中给出了计算直线与二维固定网格交线长度的算法,我们将这一算法进行推广,给出了计算直线与三维固定网格的交线长度的算法。此外,我们也给出了直线与三维固定网格的交线索引的公式。最后我们通过MATLAB实现这一算法 [10] [11] ,并用两组数值实验分别验证了算法的正确性和有效性。

2. 计算直线与网格的交线长度及索引位置

2.1. 空间坐标系的建立及图像的离散化

为后续描述方便,我们首先定义物体的参数(见表1)。

1) 假设待测物体处于一个长、宽、高分别为的矩形区域内。以矩形的中心为原点建立空间直角坐标系,对矩形区域进行离散化后得每个格点的长、宽、高分别为

Figure 1. Discretization of the object to be imaged

图1. 待测物体离散化

Table 1. Parameters of the object to be imaged

表1. 待测物体参数

2) 我们将网格中格点依轴依次从下往上逐层标号,以最底层格点为例对网格进行标号,从数字1开始,按图2的方法,直到标完对应的格点数; 按同样的方法在矩形区域内向上一层对格点进行标号,依次向上标,直到标完对应的格点数,以此完成对矩形区域的离散化。

2.2. 第i根射线所在位置描述

假设第根射线是由第个光源及探测器所确定,则该直线所在的空间直线参数方程可表示为(为参数)

Figure 2. Label of grid on xoy plane

图2. xoy平面网格标号

其中()为直线的方向向量,可以根据确定的光源和探测器位置的坐标确定。为后续描述方便,需保证按序大于0,若出现以下情况之一:(1),(2),(3),则令即可。

2.3. 直线与网格交线的长度与索引公式

1) 为确定直线与网格的交线长度,我们首先给出直线与网格的交点。设光源探测器直线与平面的交点分别为 ()。根据条件

删除越界的点,得到直线与轴网格线相交组的坐标,仍然记为。同理,可获得直线与轴网格线相交的组有效坐标值,为,与轴网格线相交的组有效坐标值,为

将三组数据按照坐标从小到大进行重排,若,则按从小到大重排,若,则按从小到大重排。得到一系列不重复交点,记为。则该直线与网格的交线长度可表示为

(1)

2) 为计算直线与网格交线的索引,可以根据交线起点找到该格点中坐标最小的顶点,并由此顶点推断索引值。具体的实现中,首先令

· 若,且,则修正;

· 若,且,则修正.

图3为几种特殊的相交情况,以此说明上述修正。

最后令

(a) (b) (c) (d)

Figure 3. Special cases of intersection of a line and fixed grid. Solid circle: starting point; Rhombus: end point; Hollow circle: point before correction; Star: corrected point. (a); (b) and; (c) and; (d) and, and

图3. 几种特殊的直线与网格相交情况。实心圆:起点;菱形:终点;空心圆:修正前的点;星形:修正后的点。(a);(b);(c);(d)

根据修正后的,可得交线所在格点的索引公式为

(2)

3. 数值实验

为验证算法的有效性和正确性,我们进行如下两组数值实验。算法在一台Dell电脑上完成,操作系统为Windows 7旗舰版(64位),处理器:Intel(R) Core(TM) i7-4790,CPU@3.60 GHz,RAM:8 GB。采用Matlab R_2013a版本。

3.1. 算法正确性验证

。分别取如下4组光源和探测器。经程序计算得到结果如表2所示。此外,我们对不同光源和探测器情况作了手动计算,经过比较,我们确定两者一致,算法的正确性得到验证。

Table 2. Correctness verification for the numerical method

表2. 程序正确性检测结果

3.2. 算法运行效率

我们模拟现实的光源和探测器分布情况,假设待测物体参数设置为

光源和探测平面同时绕待测物体螺旋上升。具体参数见表3

初始状态下光源位置设为。与之对应的探测器平面如图4所示。

光源和探测器平面同时绕轴旋转,第个光源转过角度为

光源位置为

与之相对应的探测器上各点坐标为

由于系数矩阵是一个大型稀疏矩阵,在实际的迭代过程中,如果将存入内存,是一个耗空间的做法;如果每次迭代都从外设读取数据,也是一个耗时的做法。因此我们考虑在每步迭代中,都计算一次系数矩阵。模拟实际情况,我们假设共采集次数据,形成组交线长度和索引标号。

Figure 4. Label of detector

图4. 平面探测器编号

Table 3. Parameters for sources and detectors

表3. 光源、探测器参数

采用上述程序,在MATLAB下运行,实际代码执行的时间为597.7133秒。由此,我们也可以大体估计每次迭代所用时间。根据数值实验的结果,我们发现运行效率可以接受。

4. 总结

本文以X-ray CT为例,研究重建问题中系数矩阵的生成方式,将系数矩阵元素和直线与三维固定网格相交联系起来,研究了直线与三维固定网格的交线长度与索引。最后用MATLAB程序语言实现了这一算法,并将算法应用到三维重建的系数矩阵生成上,完成一次矩阵的计算时间在10分钟以内。运算效率是可行的。

基金项目

国家自然科学基金青年基金项目(11401372)。

文章引用

张晓媛,王彩芳,殷绪导,徐兆亮. 直线与三维固定网格的交线长度与索引公式计算
Formulas of the Length and Index for the Intersection of a Line and a Fixed Grid in Three-Dimensional Space[J]. 应用数学进展, 2017, 06(04): 496-503. http://dx.doi.org/10.12677/AAM.2017.64059

参考文献 (References)

  1. 1. 冈萨雷斯. 数字图像处理[M]. 北京: 北京电子工业出版社, 2009.

  2. 2. 曾晖, 孙腊珍, 汪晓莲. CT计算机断层扫描成像实验[J]. 物理实验, 2008, 28(12): 9-12.

  3. 3. 庄天戈. CT原理与算法[M]. 上海: 上海交通大学出版社, 1992.

  4. 4. 潘晋孝. X射线CT迭代算法研究[M]. 北京: 北京大学数学科学学院, 2006.

  5. 5. 张顺利, 李卫斌, 唐高峰. 滤波反投影图像重建算法研究[J]. 咸阳师范学院学报, 2008, 23(4): 47-49.

  6. 6. 洪贤勇, 乔志伟. 用反投影滤波算法实现CT图像的ROI重建[J]. 电视技术, 2014, 38(7): 29-32.

  7. 7. 陈洪磊, 贺建峰, 刘俊卿, 马磊. 迭代图像重建中系统矩阵与重建图像质量关系研究[J]. 计算机应用, 2013, 33(1): 53-56, 68.

  8. 8. 王亮, 寿永熙, 秦俊平. 图像重建迭代算法的研究[J]. 黑龙江科技信息, 2007(21): 72-72.

  9. 9. 吴大瑞, 何钦铭. 一种简单的基于固定网格的空间直线索引算法[J]. 江南大学学报(自然科学版), 2005, 4(4): 394-396.

  10. 10. 张振东, 哈力旦·A. 基于MATLAB 的CT的图像三维重建的研究与实现[J]. 电子世界, 2013(3): 87-88.

  11. 11. 曾筝, 董芳华, 陈晓, 周宏, 周建中. 利用MATLAB实现CT断层图像的三维重建[J]. CT理论与应用研究, 2004, 13(2): 24-29.

期刊菜单