CN109360252B - 基于投影变换的锥束cl投影数据等效转换方法 - Google Patents
基于投影变换的锥束cl投影数据等效转换方法 Download PDFInfo
- Publication number
- CN109360252B CN109360252B CN201811069096.1A CN201811069096A CN109360252B CN 109360252 B CN109360252 B CN 109360252B CN 201811069096 A CN201811069096 A CN 201811069096A CN 109360252 B CN109360252 B CN 109360252B
- Authority
- CN
- China
- Prior art keywords
- projection
- scan
- image
- plane
- projection image
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 52
- 230000009466 transformation Effects 0.000 title claims abstract description 25
- 238000006243 chemical reaction Methods 0.000 title claims abstract description 22
- 238000002591 computed tomography Methods 0.000 claims abstract description 85
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 49
- 238000013507 mapping Methods 0.000 claims description 20
- 230000005855 radiation Effects 0.000 claims description 15
- 238000005516 engineering process Methods 0.000 description 8
- 238000010586 diagram Methods 0.000 description 7
- 238000003384 imaging method Methods 0.000 description 7
- PCTMTFRHKVHKIS-BMFZQQSSSA-N (1s,3r,4e,6e,8e,10e,12e,14e,16e,18s,19r,20r,21s,25r,27r,30r,31r,33s,35r,37s,38r)-3-[(2r,3s,4s,5s,6r)-4-amino-3,5-dihydroxy-6-methyloxan-2-yl]oxy-19,25,27,30,31,33,35,37-octahydroxy-18,20,21-trimethyl-23-oxo-22,39-dioxabicyclo[33.3.1]nonatriaconta-4,6,8,10 Chemical compound C1C=C2C[C@@H](OS(O)(=O)=O)CC[C@]2(C)[C@@H]2[C@@H]1[C@@H]1CC[C@H]([C@H](C)CCCC(C)C)[C@@]1(C)CC2.O[C@H]1[C@@H](N)[C@H](O)[C@@H](C)O[C@H]1O[C@H]1/C=C/C=C/C=C/C=C/C=C/C=C/C=C/[C@H](C)[C@@H](O)[C@@H](C)[C@H](C)OC(=O)C[C@H](O)C[C@H](O)CC[C@@H](O)[C@H](O)C[C@H](O)C[C@](O)(C[C@H](O)[C@H]2C(O)=O)O[C@H]2C1 PCTMTFRHKVHKIS-BMFZQQSSSA-N 0.000 description 6
- 238000011426 transformation method Methods 0.000 description 6
- 238000004364 calculation method Methods 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 230000006870 function Effects 0.000 description 3
- 230000007246 mechanism Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 239000013598 vector Substances 0.000 description 3
- 230000007423 decrease Effects 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 230000005469 synchrotron radiation Effects 0.000 description 2
- 235000014698 Brassica juncea var multisecta Nutrition 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 244000275904 brauner Senf Species 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000007429 general method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
- G06T11/008—Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明实施例提供一种基于投影变换的锥束CL投影数据等效转换方法，所述方法包括：将CL扫描中的第一投影图像转换为CT扫描中的第二投影图像；将所述CL扫描中的第一参数转换为所述CT扫描中的第二参数；根据所述第二投影图像和所述第二参数，基于CT图像重建算法进行图像重建。本发明实施例通过根据CL扫描中点源锥束X射线的传播强度特征和透视投影特性，结合CL扫描的几何关系，将CL扫描中的投影图像等效转换为CT扫描中的投影图像，将CL扫描中的几何参数等效转换为CT扫描中的几何参数，从而直接使用任何CT重建算法对CL扫描数据进行重建，方法简单，通用性强，重建精度高。
Description
技术领域
本发明实施例属于图像重建技术领域，更具体地，涉及一种基于投影变换的锥束CL投影数据等效转换方法。
背景技术
随着探测器技术与计算机技术的发展，传统CL扫描的胶片被替换为平板探测器，从而可以在扫描过程中获得一系列连续的投影数据。对于获得的数据，在早期主要使用称为Tomosyntheis的技术进行断层重建。现代CL(Computed Laminography，计算机分层扫描)技术吸收了CT(Computed Tomography，计算机断层扫描)技术，以解析或迭代的方法进行重建。相比传统方法中一次扫描只能获得一个模糊的断层，现代CL的成像速度和图像质量都有明显提升。
现代CL技术因光源几何性质的不同划分为平行束系统和锥束系统，平行束主要使用同步辐射光源，平行束CL称为SRCL(Synchrotron Radiation Computed Laminography，同步辐射计算机分层扫描)。由于锥形束是工业及科研中相对容易获得的光源，在平行束CL成功应用后，对锥束CL的研究就成为了热点。在锥束CL技术的研究过程中，出现了多种扫描轨迹，不同的扫描轨迹不但影响机械机构的设计，其重建方法也有所不同。目前所研究的轨迹大致可以分为四类：平面轨迹，有限角摆动，圆轨迹，空间轨迹。更复杂的轨迹及重建算法也仍在研究中。
2003年，M.Yang等人提出了一种倾斜回转扫描形式的CL成像技术，如图1所示，图1中圆圈表示射线源，平行四边形表示投影平面，倾斜的圆柱形表示旋转轨迹，虚线箭头的方向表示轨迹的运动方向，
表示旋转轴矢量，∑CL表示CL扫描的投影平面。该技术可以认为是任意旋转轴下的CT重建技术。传统的CT技术中，扫描回转轴平行于探测器平面，而该技术则使用与探测器有一定夹角的回转轴。对于这种轨迹，M.Yang等人提出了一种滤波反投影型的解析重建算法。这种扫描形式的主要优势是机构易于实现，且拥有极快的重建速度。圆轨迹倾斜扫描获得的投影数据信息量较多，重建质量较好，机构也相对易于实现。
在现代断层重建算法研究中，一般将CL重建问题视为特殊轨迹的CT重建问题，所以CL重建算法与现代CT领域中的各种技术密切相关。现代CL重建从CT滤波反投影(Filtered Backprojection,FBP)算法中借鉴了许多技术，这一类算法属于解析型重建算法，重建速度快，是工程中最为常用的重建手段。在锥束CT应用中，最著名的算法是基于圆轨迹扫描的近似FBP算法，一般称为FDK算法。在CL重建领域，Myagotin论述了平行束CL的FBP重建方法，Gao等人在线性扫描轨迹的基础上研究了锥束直线滤波反投影算法(LinearFBP,LFBP)；M.Yang等在锥束圆轨迹扫描中实现了一种滤波反投影重建算法。这些算法的执行效率都很高，但是重建的图像都因数据缺失而存在一定的伪影。与解析重建相比，迭代重建型算法的通用性强，重建质量好，在投影数据缺失或者噪音较大时，迭代类重建算法的重建结果明显优于滤波反投影类算法，且无论扫描轨迹多么复杂，都可以通过迭代的方式进行重建。
如果能将CT重建算法用于CL重建，一方面可以验证新技术对CL重建的作用，另一方面也可以作为一种可行的CL重建方法。然而，如果直接对CT重建算法进行修改，使之应用于CL重建，需要开发者先正确理解待修改的CT重建算法，再进行精细地编程，其中涉及了很多数学重建理论和高性能计算方面的知识，实现难度较大；由于CT重建算法很多，如果将每个CT重建算法都进行改造并用于CL重建，工作量将非常大。综上所述，亟需一种简单、通用的方法将CT重建算法用于CL重建。
发明内容
为克服上述CT重建算法用于CL重建算法需要对CT重建算法进行修改，难度大的问题或者至少部分地解决上述问题，本发明实施例提供一种基于投影变换的锥束CL投影数据等效转换方法。
根据本发明实施例的第一方面，提供一种基于投影变换的锥束CL投影数据等效转换方法，包括：
将CL扫描中的第一投影图像转换为CT扫描中的第二投影图像；
将所述CL扫描中的第一参数转换为所述CT扫描中的第二参数；
根据所述第二投影图像和所述第二参数，基于CT图像重建算法进行图像重建。
根据本发明实施例的第二个方面，还提供一种电子设备，包括：
至少一个处理器；以及
与所述处理器通信连接的至少一个存储器，其中：
所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行第一方面的各种可能的实现方式中任一种可能的实现方式所提供的基于投影变换的锥束CL投影数据等效转换方法。
本发明实施例提供一种基于投影变换的锥束CL投影数据等效转换方法，该方法通过根据CL扫描中点源锥束X射线的传播强度特征和透视投影特性，结合CL扫描的几何关系，将CL扫描中的投影图像等效转换为CT扫描中的投影图像，将CL扫描中的几何参数等效转换为CT扫描中的几何参数，从而直接使用任何CT重建算法对CL扫描数据进行重建，方法简单，通用性强，重建精度高。
附图说明
图1为现有技术中倾斜圆轨迹CL扫描示意图；
图2为本发明实施例提供的基于投影变换的锥束CL投影数据等效转换方法整体流程示意图；
图3为本发明实施例提供的基于投影变换的锥束CL投影数据等效转换方法中倾斜圆轨迹CL扫描的正视图；
图4为本发明实施例提供的基于投影变换的锥束CL投影数据等效转换方法中CL投影平面转换为CT投影平面的示意图；
图5为本发明实施例提供的基于投影变换的锥束CL投影数据等效转换方法中透视引起的几何轮廓变形示意图；
图6本发明实施例提供的基于投影变换的锥束CL投影数据等效转换方法中第二投影图像求解示意图；
图7为本发明实施例提供的基于投影变换的锥束CL投影数据等效转换方法中旋转轴的旋转角度小于0时，三维空间顺时针旋转示意图；
图8为本发明实施例提供的基于投影变换的锥束CL投影数据等效转换方法中旋转轴的旋转角度大于0时，三维空间逆时针旋转示意图；
图9为本发明实施例提供的基于投影变换的锥束CL投影数据等效转换装置整体结构示意图；
图10为本发明实施例提供的电子设备整体结构示意图。
具体实施方式
下面结合附图和实施例，对本发明实施例的具体实施方式作进一步详细描述。以下实施例用于说明本发明实施例，但不用来限制本发明实施例的范围。
在本发明实施例的一个实施例中提供一种基于投影变换的锥束CL投影数据等效转换方法，图2为本发明实施例提供的基于投影变换的锥束CL投影数据等效转换方法整体流程示意图，该方法包括：S201，将CL扫描中的第一投影图像转换为CT扫描中的第二投影图像；
其中，CL扫描为锥束CL扫描，锥束CL扫描中的锥束X射线一般使用X光管产生，通过高速电子流轰击金属靶获得，可以近似为点光源。在理想条件下，点源X射线在空间依球面传播。点源X光在真空中传播时，其强度会随着传播距离变远而减弱，本实施例将这种特性称为“平方减弱”特性。从探测器校正角度看，这种不同可以近似认为是一种探测器响应不一致，可以在校正过程中消除，所以在带校正的平板探测器成像中，一般认为平方减弱效应可以被忽略。在本实施例中均假设所有投影数据都是校正后的数据，从而忽略平方减弱效应。换言之，本实施例中认为点源X射线在真空传播时强度不变。
从几何的角度看，在锥束X射线成像的投影过程中，引发探测器产生响应的X射线都汇聚到射线源的焦点，所以锥束X射线成像具备透视投影的特征。透视投影是一种映射，按照数学语言可以描述为：在三维欧式空间R3中，预先选定点s和平面Σ，s称为焦点，Σ称为投影平面；W是R3的子集，对于任意点u∈W，都可以确定直线su与Σ的唯一交点u'，这个过程称为将点u关于s-Σ作透视投影，u'是u的投影点。将W内的所有点都透视投影到平面Σ后可以得到集合V，从而可以将透视映射记作Tp：
为了便于进行数学描述，定义焦点s为原点O，建立任意右手坐标系Oxyz。在此坐标系下，投影平面Σ的方程可以表示为∑：Ax+By+Cz＝D。此时，对于任意空间点p:(x0,y0,z0)∈R3，它关于原点O和投影平面Σ进行透视投影时，投影点p':(x0',y0',z0')的坐标可以通过求直线Op与平面Σ的交点获得。透视映射Tp可以按如下公式表达：
在锥束X射线成像中，对于任意一条射线束，它与图像重建目标abc三角形所在的平面相交于(x1,y1,z1)，在该处的强度值为q1，与Σ相交于(x2,y2,z2)，在该处的强度值为q2。如果将强度值视作图像的灰度值，那么此时就能形成两个三维空间内的平面图像I1:{(x1,y1,z1,q1)}和I2:{(x2,y2,z2,q2)}。在本发明中，假定同一束射束在真空中传播时可以认为强度不变，即q1＝q2。若取投影平面的方程为Ax+By+Cz＝D，则将q1＝q2和公式(2)组合，可以得到以下映射：
通过该映射，便可在已知abc平面内的图像I1:{(x1,y1,z1,q1)}时，直接计算出透视投影后在Σ平面内的投影图像I2＝Txp(I1)。将图像重建目标在CL扫描中的投影图像作为第一投影图像。将CL扫描下的第一投影图像进行转换为CT扫描下的第二投影图像。其中，第二投影图像根据CT扫描中的几何参数对第一投影图像进行转换获取。本实施例中CL扫描和CT扫描中的射线源相同，旋转轨迹相同。
S202，将所述CL扫描中的第一参数转换为所述CT扫描中的第二参数；
其中，第一参数为CL扫描中的参数，第二参数为CT扫描中的参数。由于CT图像重建算法不能直接用于CL扫描获取的第一投影图像的重建，只能用于CT扫描获取的第二投影图像的重建，且CT图像重建算法需要用到CT扫描中的参数。若想使用CT图像重建算法对第一投影图像重建，就需要将CL扫描参数转换为CT重建所需扫描参数。
S203，根据所述第二投影图像和所述第二参数，基于CT图像重建算法进行图像重建。
其中，CT图像重建算法为根据CT扫描获取的投影图像进行图像重建的算法。由于第二投影图像为CT扫描中的投影图像，第二参数为CT扫描中的参数，因此，可以直接使用CT图像重建算法进行图像重建。
本实施例根据CL扫描中点源锥束X射线的传播强度特征和透视投影特性，结合CL扫描的几何关系，将CL扫描中的投影图像等效转换为CT扫描中的投影图像，将CL扫描中的几何参数等效转换为CT扫描中的几何参数，从而直接使用任何CT重建算法对CL扫描数据进行重建，方法简单，通用性强，重建精度高。
在上述实施例的基础上，本实施例中将CL扫描中的第一投影图像转换为CT扫描中的第二投影图像的步骤具体包括：将CL扫描中的第一投影平面转换为CT扫描中的第二投影平面；将所述第一投影平面上的第一投影图像映射到所述第二投影平面上，获取所述第二投影平面上的第二投影图像。
其中，第一投影平面为CL扫描中的投影平面，第二投影平面为CT扫描中的投影平面。由于本实施例中CL扫描和CT扫描中的射线源相同，旋转轨迹相同，根据旋转轨迹的旋转轴的倾斜角度和CL扫描中射线源在第一投影平面上的第一投影点即可获取CT扫描中的第二投影平面。根据透视投影原理将第一投影平面上图像重建目标的第一投影图像映射到第二投影平面上，获取第二投影平面上的第二投影图像。
在上述各实施例的基础上，本实施例中将CL扫描中的第一投影平面转换为CT扫描中的第二投影平面的步骤具体包括：获取所述CL扫描中的射线源在所述第一投影平面上的第一投影点；将过所述第一投影点且与所述CL扫描中旋转轨迹的旋转轴平行的面作为第二投影平面；其中，所述CL扫描为锥束扫描，所述旋转轨迹为倾斜圆轨迹。
图3为图1中倾斜圆轨迹CL扫描的正视图，获取CL扫描中的射线源在第一投影平面上的第一投影点FCL，如图4所示，图4中θ为CL扫描中旋转轴的倾斜角度，过第一投影点F作与旋转轴平行的面，将该面作为第二投影平面ΣCT。图4中的FDD为射线源到第一投影点FCL的距离。CL扫描等效的CT扫描是一种圆轨迹CT扫描，此时投影平面偏离射线源较远，是一种大锥角的圆轨迹CT扫描。在进行CL扫描时，只需要基于透视映射原理将第一投影平面上的第一投影图像映射到第二投影平面上，获得第二投影平面上的第二投影图像，再计算出相应的CT几何参数，就可以通过CT重建算法进行重建。
在上述实施例的基础上，本实施例中将所述第一投影平面上的第一投影图像映射到所述第二投影平面上，获取所述第二投影平面上的第二投影图像的步骤具体包括：获取所述第一投影图像映射到所述第二投影平面上的梯形区域；获取所述梯形区域中平行的两条边，将平行的所述两条边中长度最短的边进行延长，获取矩形区域；根据生成所述第一投影图像的探测器中感光元件的尺寸，对所述矩形区域进行网格划分，计算每个所述网格对应的像素强度。
如图5所示，图5中的圆柱形为射线源，第一投影平面上的投影区域为矩形区域ABCD，如图5中的竖线填充区域。当矩形框内的第一投影图像映射到第二投影平面时，会发生变形，映射后的获取到第二投影平面上的梯形区域A'B'C'D'，如图5中的横线填充区域。由于A'B'平行于C'D'，所以选择延长A'B'和C'D'这两条边中较短的一个，最终形成一个包含A'B'C'D'的矩形区域PQRS，如图6所示。此时，在矩形区域PQRS内根据生成所述第一投影图像的探测器中感光元件的尺寸Psize划分网格，使每个网格的尺寸与感光元件的尺寸相同，如图6所示。计算PQRS内每个网格处对应的透视映射值，即像素强度，从而得到第二投影图像Img'。第二投影图像Img'的求解问题描述为：已知ABCD内的第一投影图像和成像几何参数，求解PQRS内的第二投影图像。
在上述实施例的基础上，本实施例中通过以下公式获取所述第一投影图像映射到所述第二投影平面上的梯形区域：
其中，(xw'，yw'，zw')为所述第一投影图像中的任一像素投影到所述第二投影平面上的第二投影点的三维坐标，(xw，yw，zw)为所述第一投影图像中任一像素的三维坐标，FDD为所述射线源和所述第一投影点之间的第一距离，θ为所述CL扫描中旋转轴的倾斜角度。
为了便于计算和分析，本实施例中将射线源作为原点O，将射线源向第一投影平面做垂直投影，得到第一投影点FCL，定义y轴平行于线段OFCL。x轴和z轴分别平行于第一投影平面相邻的两条边，且与y轴形成右手坐标系。限定CL旋转轴矢量
平行于Oyz平面，且
有正的y分量。定义
和Oxy平面的夹角为θ，在沿着-x轴进行观测时，取逆时针为正向，则有
若
则CL扫描退化为普通的CT扫描。在进行CL扫描时，若第一投影平面相对射线源的位置产生了变化，所建立的坐标系也会相应的变化。
对于第二投影平面ΣCT，其法向量为(0,-sinθ,cosθ)，且过点(0,FDDCL,0)，易得其方程为：
ΣCT:sinθy-cosθz-FDDsinθ＝0 (4)
根据透视投影公式(2)和ΣCT的方程(4)，获取将ABCD内的点(xw,yw,zw)投影到ΣCT平面上的映射公式：
根据公式(5)获取第一投影图像中各像素投影到第二投影平面上的第二投影点，将第一投影图像中各像素投影到所述第二投影平面上的投影点作为第二投影点。第二投影点组成的区域为梯形区域。
在上述各实施例的基础上，本实施例中获取所述第一投影图像映射到所述第二投影平面上的梯形区域的步骤之前还包括通过以下公式获取所述第一投影图像中任一像素的三维坐标(xw，yw，zw)：
xw＝(i-Fi)*Psize；
yw＝FDD；
zw＝(Fj-j)*Psize；
其中，i∈[1,width]，j∈[1,height]，width为所述第一投影图像的宽度，height为所述第一投影图像的高度，i为所述任一像素所在的列数，j为所述任一像素所在的行数，Fi为所述第一投影点所在的列数，Fj为所述第一投影点所在的行数，Psize为所述感光元件的尺寸。
其中，第一投影图像为二维图像，第一投影图像中的像素与感光元件的空间位置一一对应。设第一投影图像Img的宽和高分别为width和height，探测器上的感光元件的尺寸为Psize，则Img映射到三维空间后，其坐标(xw,yw,zw,qw)集合W满足：
选取第一投影图像的左上角为原点，向右为i轴，向下为j轴，二维图像可以用(i,j,Img(i,j))表示。Img(i,j)表示坐标为(i，j)的像素的强度。射线源在第一投影平面上的投影点可以使用图像空间内的坐标描述，令其为(Fi,Fj)。
在上述各实施例的基础上，本实施例中计算每个所述网格对应的像素强度的步骤具体包括：对于坐标为(xt、yt、zt)的任一所述网格，基于双线性插值算法获取所述第一投影图像中坐标
处的像素强度，将所述像素强度作为该网格对应的像素强度；其中，FDD为所述射线源和所述第一投影点之间的第一距离。
具体地，将ABCD四个顶点的坐标A(1,1),B(1,height),C(width,1),D(width,height)先代入公式(6)，得到其在三维空间内的坐标，再将得到的三维坐标代入公式(5)中就可以得到A'B'C'D'四个投影点的空间坐标。PQRS区域内像素的x坐标范围为[xP,xQ]，从数值上看，该范围和[min(xA',xD'),max(xB',xC')]相同。PQRS区域内像素的y坐标范围为[yP,yS]，从数值上看，该范围和[yA',yD']相同，令PQRS面上网格的宽度为Psize。可通过以下公式求得PQRS区域内网格的三维坐标(xt,yt,zt)中(xt,yt)的范围：
其中linspace(a,b,n)是一个函数，该函数返回一个离散的区间，该区间内的点是在连续区间[a,b]内均匀取n个点形成的(含a,b在内)。“×”运算是区间的直积，两个离散区间的直积将形成一个离散的网格点集合。将(xt,yt)代入公式(4)中平面ΣCT的方程后，计算出zt的值，从而得到PQRS内网格点的三维空间坐标集合：
T:{(xt,yt,zt)} (8)
通过插值方法获取PQRS区域中各点处对应的强度。对于PQRS区域内的点(xt,yt,zt)，通过透视映射公式(5)知，该点处的强度值和ABCD区域内的点
处相同，而
处的强度值可以通过在ABCD平面内进行二维插值得到。此时的问题相当于：已知若干点的坐标值(xw,zw)和对应的函数值f(xw,zw)＝qw，求解函数值
本实施例使用双线性插值算法进行计算。双线性插值中，在已知图像I(i,j)求解
时，
的值由距离(x0,y0)最近的四个像素值加权求和得到。在双线性插值中，认为图像的灰度在x,y方向上都是线性变化的，基于此，就很容易求出四个像素点的对应权值，进而求出插值公式：
其中，(xi,yi)为ABCD中距离(x0,y0)最近的四个点。在插值运算完成后，只需要将PQRS内每个网格点对应的值依矩阵的形式进行保存，即可得到第二投影图像Img'。
在上述各实施例的基础上，本实施例中所述第一参数包括所述CL扫描中的射线源在所述第一投影平面上第一投影点的坐标，以及所述射线源和所述第一投影点之间的第一距离；所述第二参数包括所述射线源在所述第二投影平面上第三投影点的坐标，以及所述射线源和所述第三投影点之间的第二距离。
在上述各实施例的基础上，本实施例中将所述CL扫描中的第一参数转换为所述CT扫描中的第二参数的步骤具体包括：根据所述第一距离和所述CL扫描中旋转轴的倾斜角度，获取所述第二距离；根据所述矩形区域的顶点坐标和所述第二距离，获取所述第三投影点的坐标。
CL扫描的几何参数第一距离FDD和第一投影点FCL:(Fi,Fj)按图4中的描述，CT扫描几何参数第二距离FDDCT和第三投影点FCT:(FCTi,FCTj)按图5中的描述。当θ＜0时，由于OFCT垂直于ΣCT平面，所以可以得到OFCT与y轴的夹角
得出：
FDDCT＝OFCT＝OFcosα＝FDDcosα (10)
对于(FCTi,FCTj)的计算，在第二投影平面中P点对应第二投影图像的左上角，所以(FCTi,FCTj)可以通过在第二投影平面内求FCT相对于P点的位置来得到。为了便于计算，将三维空间绕x轴顺时针旋转α，得到图7。旋转后P点到达P′：(xp′，yp′，zp′)，FCT'则恰好位于y轴的(0,FDDcosα,0)处。通过以下公式对P'的坐标进行计算：
由于三维旋转不改变点的相对关系，因此，只需在旋转后的空间内计算FCT'相对于P'的位置，就能求解出(FCTi,FCTj)，由于P'Q'R'S'平行于Oxz面，可以作为一个二维平面分析。在已知像素尺寸即感光元件尺寸为psize时，可以得出FCT'在以P'为原点的图像坐标系中的坐标(FCTi,FCTj)满足：
当θ＞0时，推理过程和θ＜0时基本相同。首先取OFCT与y轴的夹角
将三维空间绕Ox轴逆时针旋转α，如图8所示。经过同样的推理，可以得到：
此时，通过以下公式对P'的坐标进行计算：
在本发明实施例的另一个实施例中提供一种基于投影变换的锥束CL投影数据等效转换装置，该装置用于实现前述各实施例中的方法。因此，在前述基于投影变换的锥束CL投影数据等效转换方法的各实施例中的描述和定义，可以用于本发明实施例中各个执行模块的理解。图9为本发明实施例提供的基于投影变换的锥束CL投影数据等效转换装置整体结构示意图，该装置包括第一转换模块901、第二转换模块902和图像重建模块903；其中：
第一转换模块901用于将CL扫描中的第一投影图像转换为CT扫描中的第二投影图像；第二转换模块902用于将所述CL扫描中的第一参数转换为所述CT扫描中的第二参数；图像重建模块903用于根据所述第二投影图像和所述第二参数，基于CT图像重建算法进行图像重建。
在上述实施例的基础上，本实施例中第一转换模块包括第一转换子模块和映射子模块，其中：
第一转换子模块用于将CL扫描中的第一投影平面转换为CT扫描中的第二投影平面；映射子模块用于将所述第一投影平面上的第一投影图像映射到所述第二投影平面上，获取所述第二投影平面上的第二投影图像。
在上述实施例的基础上，本实施例中第一转换子模块具体用于：获取所述CL扫描中的射线源在所述第一投影平面上的第一投影点；将过所述第一投影点且与所述CL扫描中旋转轨迹的旋转轴平行的面作为第二投影平面；其中，所述CL扫描为锥束扫描，所述旋转轨迹为倾斜圆轨迹。
在上述实施例的基础上，本实施例中映射子模块具体用于：获取所述第一投影图像映射到所述第二投影平面上的梯形区域；获取所述梯形区域中平行的两条边，将平行的所述两条边中长度最短的边进行延长，获取矩形区域；根据生成所述第一投影图像的探测器中感光元件的尺寸，对所述矩形区域进行网格划分，计算每个所述网格对应的像素强度。
在上述各实施例的基础上，本实施例中映射子模块进一步用于通过以下公式获取所述第一投影图像映射到所述第二投影平面上的梯形区域：
其中，(xw'，yw'，zw')为所述第一投影图像中的任一像素投影到所述第二投影平面上的第二投影点的三维坐标，(xw，yw，zw)为所述第一投影图像中任一像素的三维坐标，FDD为所述射线源和所述第一投影点之间的第一距离，θ为所述CL扫描中旋转轴的倾斜角度。
在上述各实施例的基础上，本实施例中映射子模块还用于通过以下公式获取所述第一投影图像中任一像素的三维坐标(xw，yw，zw)：
xw＝(i-Fi)*Psize；
yw＝FDD；
zw＝(Fj-j)*Psize；
其中，i∈[1,width]，j∈[1,height]，width为所述第一投影图像的宽度，height为所述第一投影图像的高度，i为所述任一像素所在的列数，j为所述任一像素所在的行数，Fi为所述第一投影点所在的列数，Fj为所述第一投影点所在的行数，Psize为所述感光元件的尺寸。
在上述实施例的基础上，本实施例中映射子模块进一步用于：对于坐标为(xt、yt、zt)的任一所述网格，基于双线性插值算法获取所述第一投影图像中坐标
处的像素强度，将所述像素强度作为该网格对应的像素强度；其中，FDD为所述射线源和所述第一投影点之间的第一距离。
在上述各实施例的基础上，本实施例中所述第一参数包括所述CL扫描中的射线源在所述第一投影平面上第一投影点的坐标，以及所述射线源和所述第一投影点之间的第一距离；所述第二参数包括所述射线源在所述第二投影平面上第三投影点的坐标，以及所述射线源和所述第三投影点之间的第二距离。
在上述各实施例的基础上，本实施例中第二转换模块具体用于：根据所述第一距离和所述CL扫描中旋转轴的倾斜角度，获取所述第二距离；根据所述矩形区域的顶点坐标和所述第二距离，获取所述第三投影点的坐标。
本实施例通过根据CL扫描中点源锥束X射线的传播强度特征和透视投影特性，结合CL扫描的几何关系，将CL扫描中的投影图像等效转换为CT扫描中的投影图像，将CL扫描中的几何参数等效转换为CT扫描中的几何参数，从而直接使用任何CT重建算法对CL扫描数据进行重建，方法简单，通用性强，重建精度高。
本实施例提供一种电子设备，图10为本发明实施例提供的电子设备整体结构示意图，该设备包括：至少一个处理器1001、至少一个存储器1002和总线1003；其中，
处理器1001和存储器1002通过总线1003完成相互间的通信；
存储器1002存储有可被处理器1001执行的程序指令，处理器调用程序指令能够执行上述各方法实施例所提供的方法，例如包括：将CL扫描中的第一投影图像转换为CT扫描中的第二投影图像；将所述CL扫描中的第一参数转换为所述CT扫描中的第二参数；根据所述第二投影图像和所述第二参数，基于CT图像重建算法进行图像重建。
本实施例提供一种非暂态计算机可读存储介质，非暂态计算机可读存储介质存储计算机指令，计算机指令使计算机执行上述各方法实施例所提供的方法，例如包括：将CL扫描中的第一投影图像转换为CT扫描中的第二投影图像；将所述CL扫描中的第一参数转换为所述CT扫描中的第二参数；根据所述第二投影图像和所述第二参数，基于CT图像重建算法进行图像重建。
本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。
以上所描述的电子设备实施例仅仅是示意性的，其中作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性的劳动的情况下，即可以理解并实施。
通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到各实施方式可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件。基于这样的理解，上述技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在计算机可读存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行各个实施例或者实施例的某些部分的方法。
最后，本申请的方法仅为较佳的实施方案，并非用于限定本发明实施例的保护范围。凡在本发明实施例的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明实施例的保护范围之内。
Claims (7)
1.一种基于投影变换的锥束CL投影数据等效转换方法，其特征在于，包括：
将CL扫描中的第一投影图像转换为CT扫描中的第二投影图像；
将所述CL扫描中的第一参数转换为所述CT扫描中的第二参数；
根据所述第二投影图像和所述第二参数，基于CT图像重建算法进行图像重建；
其中，将CL扫描中的第一投影图像转换为CT扫描中的第二投影图像的步骤具体包括：
将CL扫描中的第一投影平面转换为CT扫描中的第二投影平面；
将所述第一投影平面上的第一投影图像映射到所述第二投影平面上，获取所述第二投影平面上的第二投影图像；
其中，将CL扫描中的第一投影平面转换为CT扫描中的第二投影平面的步骤具体包括：
获取所述CL扫描中的射线源在所述第一投影平面上的第一投影点；
将过所述第一投影点且与所述CL扫描中旋转轨迹的旋转轴平行的面作为第二投影平面；其中，所述CL扫描为锥束扫描，所述旋转轨迹为倾斜圆轨迹；
其中，将所述第一投影平面上的第一投影图像映射到所述第二投影平面上，获取所述第二投影平面上的第二投影图像的步骤具体包括：
获取所述第一投影图像映射到所述第二投影平面上的梯形区域；
获取所述梯形区域中平行的两条边，将平行的所述两条边中长度最短的边进行延长，获取矩形区域；
根据生成所述第一投影图像的探测器中感光元件的尺寸，对所述矩形区域进行网格划分，计算每个所述网格对应的像素强度。
2.根据权利要求1所述的方法，其特征在于，通过以下公式获取所述第一投影图像映射到所述第二投影平面上的梯形区域：
其中，(xw'，yw'，zw')为所述第一投影图像中的任一像素投影到所述第二投影平面上的第二投影点的三维坐标，(xw，yw，zw)为所述第一投影图像中任一像素的三维坐标，FDD为射线源和所述第一投影点之间的第一距离，θ为所述CL扫描中旋转轴的倾斜角度。
3.根据权利要求2所述的方法，其特征在于，获取所述第一投影图像映射到所述第二投影平面上的梯形区域的步骤之前还包括通过以下公式获取所述第一投影图像中任一像素的三维坐标(xw，yw，zw)：
xw＝(i-Fi)*Psize；
yw＝FDD；
zw＝(Fj-j)*Psize；
其中，i∈[1,width]，j∈[1,height]，width为所述第一投影图像的宽度，height为所述第一投影图像的高度，i为所述任一像素所在的列数，j为所述任一像素所在的行数，Fi为所述第一投影点所在的列数，Fj为所述第一投影点所在的行数，Psize为所述感光元件的尺寸。
4.根据权利要求1所述的方法，其特征在于，计算每个所述网格对应的像素强度的步骤具体包括：
对于坐标为(xt、yt、zt)的任一所述网格，基于双线性插值算法获取所述第一投影图像中坐标
处的像素强度，将所述像素强度作为该网格对应的像素强度；其中，FDD为射线源和所述第一投影点之间的第一距离。
5.根据权利要求1所述的方法，其特征在于，所述第一参数包括所述CL扫描中的射线源在所述第一投影平面上第一投影点的坐标，以及所述射线源和所述第一投影点之间的第一距离；
所述第二参数包括所述射线源在所述第二投影平面上第三投影点的坐标，以及所述射线源和所述第三投影点之间的第二距离。
6.根据权利要求5所述的方法，其特征在于，将所述CL扫描中的第一参数转换为所述CT扫描中的第二参数的步骤具体包括：
根据所述第一距离和所述CL扫描中旋转轴的倾斜角度，获取所述第二距离；
根据所述矩形区域的顶点坐标和所述第二距离，获取所述第三投影点的坐标。
7.一种电子设备，其特征在于，包括：
至少一个处理器、至少一个存储器和总线；其中，
所述处理器和存储器通过所述总线完成相互间的通信；
所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行如权利要求1至6任一所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811069096.1A CN109360252B (zh) | 2018-09-13 | 2018-09-13 | 基于投影变换的锥束cl投影数据等效转换方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811069096.1A CN109360252B (zh) | 2018-09-13 | 2018-09-13 | 基于投影变换的锥束cl投影数据等效转换方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109360252A CN109360252A (zh) | 2019-02-19 |
| CN109360252B true CN109360252B (zh) | 2020-08-14 |
Family
ID=65350623
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811069096.1A Active CN109360252B (zh) | 2018-09-13 | 2018-09-13 | 基于投影变换的锥束cl投影数据等效转换方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109360252B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113052929B (zh) * | 2021-03-09 | 2024-07-12 | 重庆大学 | 一种基于投影视角加权的直线扫描cl重建方法 |
| CN113409412B (zh) * | 2021-06-02 | 2023-01-10 | 北京航空航天大学 | 偏置扫描模式下的大视场cl重建方法和装置、设备及介质 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102123664A (zh) * | 2008-08-13 | 2011-07-13 | 皇家飞利浦电子股份有限公司 | 使用基于校准体模的旋转中心建立算法在不理想等中心3d旋转x射线扫描器系统中进行环形伪影校正的校准方法 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105596022B (zh) * | 2006-02-27 | 2020-05-22 | 罗切斯特大学 | 锥束ct动态成像的方法和设备 |
| EP2024935A1 (en) * | 2006-05-11 | 2009-02-18 | Philips Intellectual Property & Standards GmbH | Method and apparatus for reconstructing an image |
| US8938111B2 (en) * | 2010-01-13 | 2015-01-20 | Fei Company | Computed tomography imaging process and system |
| CN102004111B (zh) * | 2010-09-28 | 2012-09-19 | 北京航空航天大学 | 一种倾斜多锥束直线轨迹ct成像方法 |
| CN107796835B (zh) * | 2017-10-20 | 2021-05-25 | 北京航空航天大学 | 一种x射线柱面三维锥束计算机层析成像方法及装置 |
-
2018
- 2018-09-13 CN CN201811069096.1A patent/CN109360252B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102123664A (zh) * | 2008-08-13 | 2011-07-13 | 皇家飞利浦电子股份有限公司 | 使用基于校准体模的旋转中心建立算法在不理想等中心3d旋转x射线扫描器系统中进行环形伪影校正的校准方法 |
Non-Patent Citations (4)
| Title |
|---|
| Automatic measurement of rotation center for laminography scanning system without dedicated phantoms;Min Yang,etc;《Journal of Electronic Imaging》;20140930;第23卷(第5期);第1页至第12页 * |
| 一种基于投影变换的旋转型CL重建算法;王敬雨等;《CT理论与应用研究》;20171031;第26卷(第5期);第575页至582页 * |
| 基于非对称双边截断投影数据的CT重建方法;韩旭等;《全国射线数字成像与CT新技术》;20160930;第1页至第7页 * |
| 锥束射线倾斜扫描三维层析算法;杨民等;《上海交通大学学报》;20060731;第40卷(第7期);第1079页至1083页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109360252A (zh) | 2019-02-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP4409043B2 (ja) | トモシンセシスｘ線イメージング・システムにより取得される画像データを再構成するコンピュータ・プログラム及び装置 | |
| Thibault et al. | A three‐dimensional statistical approach to improved image quality for multislice helical CT | |
| JP4695867B2 (ja) | 非一様ビュー加重式トモシンセシス方法及び装置 | |
| US9861332B2 (en) | Tomographic image generation device and method, and recording medium | |
| Karolczak et al. | Implementation of a cone‐beam reconstruction algorithm for the single‐circle source orbit with embedded misalignment correction using homogeneous coordinates | |
| JP4342164B2 (ja) | コンピュータ断層撮影装置 | |
| CN110057847B (zh) | Tr层析扫描投影重排方法及装置 | |
| KR101591381B1 (ko) | Ct 촬영에서의 금속에 의한 잡음 및 오류 감쇄방법 | |
| CN109360252B (zh) | 基于投影变换的锥束cl投影数据等效转换方法 | |
| AU2017203626A1 (en) | A method and apparatus for motion correction in CT imaging | |
| Niebler et al. | Projection‐based improvement of 3D reconstructions from motion‐impaired dental cone beam CT data | |
| JP6881827B2 (ja) | Ｘ線撮像におけるグリッドアーチファクトの学習に基づく補正 | |
| Sunnegårdh et al. | Regularized iterative weighted filtered backprojection for helical cone‐beam CT | |
| Renders et al. | Adjoint image warping using multivariate splines with application to four‐dimensional computed tomography | |
| Li et al. | Motion correction for robot-based x-ray photon-counting CT at ultrahigh resolution | |
| US20170365076A1 (en) | Virtual projection image method | |
| Preuhs et al. | Fast epipolar consistency without the need for pseudo matrix inverses | |
| Bruder et al. | Compensation of skull motion and breathing motion in CT using data-based and image-based metrics, respectively | |
| CA2825704A1 (en) | Optimized implementation of back projection for computed tomography (ct) | |
| Thies et al. | A gradient-based approach to fast and accurate head motion compensation in cone-beam CT | |
| JP6243219B2 (ja) | 画像生成装置および放射線断層撮影装置並びにプログラム | |
| Chen et al. | General rigid motion correction for computed tomography imaging based on locally linear embedding | |
| Jun et al. | Alignment theory of parallel-beam computed tomography image reconstruction for elastic-type objects using virtual focusing method | |
| JP7178621B2 (ja) | 画像処理装置及び画像処理方法 | |
| WO2020209312A1 (ja) | 検査装置及び検査方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |