专利名称:Ct投影数据三维解析模拟方法
技术领域：
本发明涉及CT投影图像重建的方法，特别涉及CT投影图像重建的模拟方法。
背景技术：
CT技术广泛应用于工业无损检测(NDT)、医学诊断和科学研究领域。CT重建算法的研究过程中必须使用CT投影数据。
CT投影数据的来源可以有两种1)从实际CT系统采集；2)通过计算机模拟计算产生。
由于实际CT系统价格昂贵，多数研究单位并没有实际的CT系统。另外，为了研究重建算法的性能需要特定模型的投影数据，在实际CT系统中采集投影数据时必须定做相应的模型，开销很大，而且受实际制造能力的制约，有些模型无法制造出来。基于上述原因，计算机模拟计算投影数据在CT重建算法研究过程中被广泛采用。
计算机模拟计算投影数据的方法有两类1)解析模拟方法；2)数值模拟方法。解析模拟方法通过数学表达式直接计算获得模拟投影数据，所需的模型定义直接写在计算机程序代码中，如果需要更改模型则必须重新编写模拟投影计算程序，但其模拟投影数据是理想的，不含有离散化噪声。数值模拟方法则将模型作为像素的集合处理，模拟过程为计算各个像素并将结果合并以获得模拟投影数据，模型可由一个个的像素定义，更改模型无需更改模拟投影计算程序，因此使用非常方便。但是，数值模拟方法有两点不足1)计算速度慢；2)投影过程中引入了离散化噪声。减小离散化噪声的方法是减小像素尺寸，但这样则会导致计算时间的成三次方倍增长。对于CT算法研究而言，模拟投影数据中不希望有离散化噪声，实际选择的模拟投影计算方法则需要根据计算时间、对离散化噪声的容忍程度、模型更改的方便性这三者之间折中，选择合适的计算方法，目前条件下解析模拟方法的计算速度通常为数值模拟方法的几十至几百倍。
发明内容(一)要解决的技术问题本发明的目的是提供一种计算速度快、便于修改模型的CT投影数据三维解析模拟方法，从而能够方便、快速地获得各种模型的无离散化噪声的模拟投影数据。
(二)技术方案为了达到上述目的，本发明采取以下方案包括1)确定模型定义方法；2)确定扫描方式定义方法；3)在计算机中设置解析模拟程序；4)由1)、2)、3)产生CT投影数据pjk。
其中，所述模型定义方法包括A通过多个指定了各自衰减系数的基本几何体的组合定义一个模型；B设一个模型由i个基本几何体Mi组成，各个基本几何体的衰减系数为mi；C通过计算机交互界面进行模型定义。
其中，所述扫描方式定义方法包括a完整的模拟投影过程由j个投影构成；b为每个投影指定此投影的射线源中心位置
探测器中心位置 c指定探测器的探测单元数k及各个探测器单元中心
相对探测器中心位置
的几何关系；d指定每个探测器单元接收的射线数量l、各条射线起点
相对射线源中心位置
的几何关系及各条射线与探测器交点
相对探测器单元中心
的几何关系、各条射线对探测器单元的模拟投影数值的贡献权重wjkl，权重应是归一化的，即
e通过计算机交互界面进行扫描方式的定义。
其中，所述解析模拟程序包括a)将设置的待模拟物理因素分为对射线位置的影响和对投影值的影响两类。根据对射线位置的影响调整射线的起点和终点坐标，根据对投影值的影响调整最终模拟投影值。
b)通过坐标转换将射线起点、终点坐标转换到几何体局部坐标系，在几何体局部坐标系中计算射线与几何体的交点坐标。
c)对直接叠加组合方式，通过几何体局部坐标系中的坐标计算交线长度，随后计算衰减系数。
d)对替换组合方式，将交点坐标转换到系统坐标系，在系统坐标系中计算两个交点到射线起点的距离，随后计算衰减系数。
其中，所述CT投影数据pjk包括A)一条射线的模拟投影数值pjkl，为线段
上所有衰减系数m的积分，即pjkl=O`rrjkldu]]>B)各个探测器单元的模拟投影数值，pjk为 C)pjkl的计算方法和模型的衰减系数有关，模型的衰减系数由各个几何体的衰减系数按一定的方式组合得到，不同的组合方式必须使用不同的计算方法计算，pjkl，其中，所述组合方式包括各个几何体衰减系数的直接累加，即设模型中的一点
位于组成此模型的m个几何体内，则模型该点的衰减系数为
这种组合方式下，pjkl的计算方法为
先计算线段
与组成模型的各个几何体的交线长度Lijkl，然后计算模拟投影数值pjkl 其中，所述组合方式包括各个几何体衰减系数的替换组合，即设模型中的一点
位于组成此模型的m个几何体内，则模型该点的衰减系数为mrr=mm,]]>这种组合方式下pjkl的计算方法为设线段
与组成模型的i个几何体中的N个几何体有两个交点，和第n个几何体的两个交点位于第n层，到
的距离分别是ajkln和bjkln，设ajkln＜bjkln。ajkl0＝0，bjkl0=|rrjkl|,]]>m0＝0。第0层为线段
，各层的线段都可能被更高层的交点分割为更短的线段，各点的衰减系数值为其所属的最高层线段对应的几何体的衰减系数值。由此可计算模拟投影数值，pjkl。
(三)有益效果1)由于能提供给用于以交互的方式改变模型定义与扫描参数定义，实现在无需重新编写模拟投影程序的情况下根据需要定义各种模型与扫描参数，从而方便快捷地获得解析模拟投影数据。2)计算速度快，能获得各种模型的无离散化噪声的模拟投影数据。
图1是本发明各个几何体衰减系数的替换组合时一条射线的模拟投影值计算流程图；图2是本发明解析模拟投影过程示意图；图3是本发明模型定义交互界面示意图；图4是本发明模拟投影参数设置交互界面示意图。
图5是本发明解析模拟投影程序计算流程图。
1)具体实施方式
以下实施例用于说明本发明，但不用来限制本发明的范围。
本发明采实施时，包括步骤1)确定模型定义方法；2)确定扫描方式定义方法；3)在计算机中设置解析模拟程序；4)由步骤1)、步骤2)、步骤3)产生CT投影数据pjk。
如图2所示，本发明的典型实施例包括计算机交互界面及解析模拟程序对用户输入的处理方法。实施例中使用两个坐标系统系统坐标系与几何体局部坐标系。模型中心位于系统坐标系原点，射线源位置、探测器位置都在系统坐标系中定义。几何体中心位于几何体局部坐标系原点，对称轴与坐标轴重合。
模型定义的交互界面如图3所示，用户确定组成模型的几何体中心在系统坐标系中的坐标(xi，yi，zi)及几何体的对称轴在系统坐标系中的倾斜角度qi、fi，同时确定几何体衰减系数mi及其组合方式。
模拟投影参数设置的交互界面如图4所示，用户确定扫描方式及要模拟的物理影响因素。
解析模拟程序的流程如图5所示，对用户输入的处理方法为1)将要模拟的物理影响因素分解为两类对射线位置的影响(如各种偏移等)及对射线上投影值的影响(如噪声等)。
2)根据扫描参数设置及影响射线位置的物理因素设置计算得到各条射线的起点坐标(xsjkl，ysjkl，zsjkl)和终点坐标(xdjkl，ydjkl，zdjkl)。
3)将起点坐标和终点坐标通过坐标转换得到各条射线在各个几何体的局部坐标系中的起点坐标
和终点坐标
4)在各个几何体的局部坐标系中计算射线与几何体的交点坐标没有两个交点(无交点或仅有一个交点)；或者有两个交点 和
如果几何体的衰减系数按直接累加方式组合，根据两点间距离公式计算射线和各几何体交线长度Lijkl，然后即可计算得到各条射线的投影值pjkl。pjkl的计算方法为先计算线段
与组成模型的各个几何体的交线长度Lijkl，然后计算模拟投影数值pjkl
如果几何体的衰减系数按替换方式组合，则将交点坐标通过坐标转换得到其在系统坐标系中的坐标(xaijkl，yaijkl，zaijkl)和(xbijkl，ybijkl，zbijkl)，再计算交点到射线起点的距离aijkl和bijkl，然后即可计算得到各条射线的投影值pjkl。
相应的计算公式为 aijkl=(xaijkl-xsjkl)2+(yaijkl-ysjkl)2+(zaijkl-zsjkl)2,]]>bijkl=(xbijkl-xsjkl)2+(ybijkl-ysjkl)2+(zbijkl-zsjkl)2.]]>各条射线对探测器单元的模拟投影数值的贡献权重wjkl可选择为wjkl=1l]]>如图1所示，上述替换方式组合下，pjkl的计算方法为设线段
与组成模型的i个几何体中的N个几何体有两个交点，和第n个几何体的两个交点位于第n层，到
的距离分别是ajkln和bjkln，设ajkln＜bjkln。ajkl0＝0，bjkl0=|rrjkl|,]]>m0＝0。第0层为线段
各层的线段都可能被更高层的交点分割为更短的线段，各点的衰减系数值为其所属的最高层线段对应的几何体的衰减系数值。由此按图1所示的方法计算模拟投影数值pjkl。将所有的交点按到
的距离从小到大的顺序排列，然后以此查找各位于最高层的线段。使用一个变量pjkl保存模拟投影值，初值置为0。使用变量Sc、Ec保存当前正在处理的线段起点、终点位置，变量nc保存当前处理的线段所在层。
Sc从0开始，按以下方法查找对应的Ec及nc如果当前处理的线段起点Sc为ajkln类点(即第n层中距离
较近的点)，则先计算Ec，Ec为所有高于当前层nc且大于Sc的ajkln类点和当前层对应的bjkln类点(即第n层中距离
较远的点)bjklnc中最小的点，然后重新计算当前层nc，将nc置为Ec所在层；如果当前处理的线段起点Sc为bjkln类点，则先重新计算当前层nc，将nc置为满足对应的ajkln类点小于Sc、对应的bjkln类点小于Ec(此时为前一次处理的线段终点)且本身层数小于nc的所有层中的最大值，然后计算当前处理的线段终点Ec，Ec为所有高于当前层nc且大于Sc的ajkln类点和当前层对应的bjkln类点(即第n层中距离
较远的点)bjklnc中最小的点。随后查找当前层nc对应的衰减系数mnc，然后将(Ec-Sc)′mnc累加到pjkl中，这样就完成了一次线段处理。随后将当前处理线段起点置为前一次处理的线段终点，然后重复以上的线段处理过程，直到起点到达所有点中的最大值bjkl0时即完成计算。
权利要求
1.CT投影数据三维解析模拟方法，其特征在于包括1)确定模型定义方法；2)确定扫描方式定义方法；3)在计算机中设置解析模拟程序；4)由1)、2)、3)产生CT投影数据pjk。
2.如权利要求
1所述的CT投影数据三维解析模拟方法，其特征在于所述模型定义方法包括A通过多个指定了各自衰减系数的基本几何体的组合定义一个模型；B设一个模型由i个基本几何体Mi组成，各个基本几何体的衰减系数为mi；C通过计算机交互界面进行模型定义。
3.如权利要求
1所述的CT投影数据三维解析模拟方法，其特征在于所述扫描方式定义方法包括a完整的模拟投影过程由j个投影构成；b为每个投影指定此投影的射线源中心位置
、探测器中心位置 c指定探测器的探测单元数k及各个探测器单元中心
相对探测器中心位置
的几何关系；d指定每个探测器单元接收的射线数量l、各条射线起点
相对射线源中心位置
的几何关系及各条射线与探测器交点
相对探测器单元中心
的几何关系、各条射线对探测器单元的模拟投影数值的贡献权重wjkl，权重应是归一化的，即
e通过计算机交互界面进行扫描方式的定义。
4.如权利要求
1所述的CT投影数据三维解析模拟方法，其特征在于所述解析模拟程序包括a)将设置的待模拟物理因素分为对射线位置的影响和对投影值的影响两类，根据对射线位置的影响调整射线的起点和终点坐标，根据对投影值的影响调整最终模拟投影值；b)通过坐标转换将射线起点、终点坐标转换到几何体局部坐标系，在几何体局部坐标系中计算射线与几何体的交点坐标；c)对直接叠加组合方式，通过几何体局部坐标系中的坐标计算交线长度，随后计算衰减系数；d)对替换组合方式，将交点坐标转换到系统坐标系，在系统坐标系中计算两个交点到射线起点的距离，随后计算衰减系数。
5.如权利要求
1所述的CT投影数据三维解析模拟方法，其特征在于所述CT投影数据pjk包括A)一条射线的模拟投影数值pjkl为线段
上所有衰减系数m的积分，即pjkl=O`rrjkldu]]>B)各个探测器单元的模拟投影数值pjk为 C)pjkl的计算方法和模型的衰减系数有关，模型的衰减系数由各个几何体的衰减系数按一定的方式组合得到，不同的组合方式必须使用不同的计算方法计算pjkl。
6.如权利要求
4所述的CT投影数据三维解析模拟方法，其特征在于所述组合方式包括各个几何体衰减系数的直接累加，即设模型中的一点
位于组成此模型的m个几何体内，则模型该点的衰减系数为
，这种组合方式下pjkl的计算方法为先计算线段
与组成模型的各个几何体的交线长度Lijkl，然后计算模拟投影数值pjkl
7.如权利要求
4所述的CT投影数据三维解析模拟方法，其特征在于所述组合方式包括各个几何体衰减系数的替换组合，即设模型中的一点
位于组成此模型的m个几何体内，则模型该点的衰减系数为mrr=mm]]>，这种组合方式下pjkl的计算方法为设线段
与组成模型的i个几何体中的N个几何体有两个交点，和第n个几何体的两个交点位于第n层，到
的距离分别是ajkln和bjkln，设ajkln＜bjkln，ajkl0＝0，bjklo=|rrjkl|]]>，m0＝0，第0层为线段
，各层的线段都可能被更高层的交点分割为更短的线段，各点的衰减系数值为其所属的最高层线段对应的几何体的衰减系数值，由此可计算模拟投影数值pjkl。
专利摘要
本发明涉及CT投影图像重建的方法。本发明公开的一种CT投影数据三维解析模拟方法，包括1)确定模型定义方法；2)确定扫描方式定义方法；3)在计算机中设置解析模拟程序；4)由1)、2)、3)产生CT投影数据p
文档编号G06T11/00GK1996391SQ200510135927
公开日2007年7月11日 申请日期2005年12月31日
发明者康克军, 张丽, 陈志强, 邢宇翔, 唐杰, 程建平, 李元景 申请人:清华大学, 清华同方威视技术股份有限公司导出引文BiBTeX, EndNote, RefMan