CN111643104B - Ct散射校正方法及系统 - Google Patents
Ct散射校正方法及系统 Download PDFInfo
- Publication number
- CN111643104B CN111643104B CN202010133896.6A CN202010133896A CN111643104B CN 111643104 B CN111643104 B CN 111643104B CN 202010133896 A CN202010133896 A CN 202010133896A CN 111643104 B CN111643104 B CN 111643104B
- Authority
- CN
- China
- Prior art keywords
- scattering
- energy spectrum
- value
- correction
- detector
- 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
- 238000012937 correction Methods 0.000 title claims abstract description 208
- 238000000034 method Methods 0.000 title claims abstract description 76
- 238000001228 spectrum Methods 0.000 claims abstract description 226
- 238000005315 distribution function Methods 0.000 claims abstract description 66
- 230000014509 gene expression Effects 0.000 claims description 83
- 238000013170 computed tomography imaging Methods 0.000 claims description 44
- 230000005540 biological transmission Effects 0.000 claims description 16
- 238000001914 filtration Methods 0.000 claims description 16
- 230000003595 spectral effect Effects 0.000 claims description 16
- 230000000694 effects Effects 0.000 claims description 14
- 238000001514 detection method Methods 0.000 claims description 8
- 230000005855 radiation Effects 0.000 claims description 8
- 230000033001 locomotion Effects 0.000 claims description 5
- 230000008859 change Effects 0.000 claims description 4
- 238000013507 mapping Methods 0.000 abstract description 20
- 230000008569 process Effects 0.000 abstract description 15
- 238000004364 calculation method Methods 0.000 abstract description 9
- 238000005259 measurement Methods 0.000 abstract description 8
- 238000002591 computed tomography Methods 0.000 description 49
- 238000004422 calculation algorithm Methods 0.000 description 12
- 238000010586 diagram Methods 0.000 description 6
- 239000000126 substance Substances 0.000 description 6
- VWDWKYIASSYTQR-UHFFFAOYSA-N sodium nitrate Chemical compound [Na+].[O-][N+]([O-])=O VWDWKYIASSYTQR-UHFFFAOYSA-N 0.000 description 5
- 238000013459 approach Methods 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 238000002083 X-ray spectrum Methods 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 2
- 230000002238 attenuated effect Effects 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- MUKYLHIZBOASDM-UHFFFAOYSA-N 2-[carbamimidoyl(methyl)amino]acetic acid 2,3,4,5,6-pentahydroxyhexanoic acid Chemical compound NC(=N)N(C)CC(O)=O.OCC(O)C(O)C(O)C(O)C(O)=O MUKYLHIZBOASDM-UHFFFAOYSA-N 0.000 description 1
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 description 1
- RYGMFSIKBFXOCR-UHFFFAOYSA-N Copper Chemical compound [Cu] RYGMFSIKBFXOCR-UHFFFAOYSA-N 0.000 description 1
- ZOKXTWBITQBERF-UHFFFAOYSA-N Molybdenum Chemical compound [Mo] ZOKXTWBITQBERF-UHFFFAOYSA-N 0.000 description 1
- 229910052774 Proactinium Inorganic materials 0.000 description 1
- BQCADISMDOOEFD-UHFFFAOYSA-N Silver Chemical compound [Ag] BQCADISMDOOEFD-UHFFFAOYSA-N 0.000 description 1
- ATJFFYVFTNAWJD-UHFFFAOYSA-N Tin Chemical compound [Sn] ATJFFYVFTNAWJD-UHFFFAOYSA-N 0.000 description 1
- 238000000333 X-ray scattering Methods 0.000 description 1
- 241000159213 Zygophyllaceae Species 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000004397 blinking Effects 0.000 description 1
- 239000002041 carbon nanotube Substances 0.000 description 1
- 229910021393 carbon nanotube Inorganic materials 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 229910052802 copper Inorganic materials 0.000 description 1
- 239000010949 copper Substances 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 229910052745 lead Inorganic materials 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 229910052750 molybdenum Inorganic materials 0.000 description 1
- 239000011733 molybdenum Substances 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 210000000056 organ Anatomy 0.000 description 1
- 238000004806 packaging method and process Methods 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 235000009165 saligot Nutrition 0.000 description 1
- 229910052709 silver Inorganic materials 0.000 description 1
- 239000004332 silver Substances 0.000 description 1
- 239000000758 substrate Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 229910052718 tin Inorganic materials 0.000 description 1
- 239000011135 tin Substances 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
- A61B6/5258—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
- A61B6/03—Computed tomography [CT]
- A61B6/032—Transmission computed tomography [CT]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/44—Constructional features of apparatus for radiation diagnosis
- A61B6/4411—Constructional features of apparatus for radiation diagnosis the apparatus being modular
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Medical Informatics (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- Biophysics (AREA)
- High Energy & Nuclear Physics (AREA)
- Veterinary Medicine (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Optics & Photonics (AREA)
- Pathology (AREA)
- Public Health (AREA)
- Physics & Mathematics (AREA)
- Heart & Thoracic Surgery (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Pulmonology (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Apparatus For Radiation Diagnosis (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
Abstract
一种CT散射校正方法及系统，该CT散射校正方法包括：基于预先建立的投影数据值与散射值映射关系查找表和测量的投影值得到散射估计值，映射关系中的所述散射值为经过能谱校正后的散射值；或者基于能谱物理对散射值的分布函数进行校正后迭代求解得到散射估计值。上述方法通过在散射估计过程中加入能谱校正的因素，基于能谱硬化对于散射分布影响的不同假设，按照迭代求解的方式或者基于散射遮挡硬件测量的方式可以得到能谱校正后的散射估计值，对于硬件的要求不高，同时还具有迭代效率高、计算准确度较高的优点。
Description
技术领域
本公开属于辐射成像技术领域，涉及一种CT散射校正方法及系统。
背景技术
计算机断层成像(CT，Computed Tomography)是使用X射线束和探测器扫描物体的某一截面(断层)，利用人体器官或组织对X射线的吸收系数不同的特性获得反映此截面物理或化学特性的投影数据集合，通过计算机运算获得截面上任意位置的参数值，并由此得到断层图像。
现有CT系统的主要部件包括：X射线源、探测器、以及旋转装置。在CT成像过程中，X射线源与探测器围绕物体进行相对旋转运动，以此获得不同旋转角度下的CT数据，也称CT投影值。在CT成像中，X射线的衰减服从指数规律，CT投影值是间接获取的，通常需要“负对数”预处理。
光子散射是CT自诞生以来就存在的影响CT图像质量的基础物理挑战，会造成图像伪影和CT值不准确等问题。提高CT成像性能的最核心问题之一是去除或减少射线散射。去除散射的方法大致可以分为：基于硬件的直接散射测量和基于算法的散射估计两大类，基于硬件的直接散射测量的方法中的硬件可以利用散射遮挡块/散射遮挡条等，基于算法的散射估计的方法中的算法包括：基于物理的解析模型玻尔兹曼输运方程/蒙特卡洛计算、投影域卷积滤波等。一般而言，直接测量类散射校正方法精度高，但对硬件有额外要求，且往往需要二次扫描，可能增加剂量；而算法估计类方法对硬件没有额外要求，也不需要二次扫描，但校正效果可能差一些，或者计算复杂度高。
因此，亟需提出一种对于硬件要求不高且计算较为简单的散射校正方法。
发明内容
(一)要解决的技术问题
本公开提供了一种CT散射校正方法及系统，以至少部分解决以上所提出的技术问题。
(二)技术方案
根据本公开的一个方面，提供了一种CT散射校正方法，包括：基于预先建立的投影数据值与散射值映射关系查找表和测量的投影值得到散射估计值，映射关系中的所述散射值为经过能谱校正后的散射值；或者基于能谱物理对散射值的分布函数进行校正后迭代求解得到散射估计值；
其中，所述基于能谱物理对散射值的分布函数进行校正后迭代求解得到散射估计值，包括：将被扫描物体放于多光源CT成像系统中，获取至少两个不同位置的光源在探测器同一像素点处的投影数据；对至少两个不同位置的光源在探测器同一像素点处的投影数据进行能谱校正，获得能谱校正后探测器像素探测得到的总光子数的表达式，该表达式中包括：经过能谱校正后散射值的分布函数，其中所述经过能谱校正后散射值的分布函数基于能谱硬化对于散射分布的影响确定；以及根据能谱校正后散射值的分布函数求解出能谱校正后散射值的表达式，并根据能谱校正后散射值的表达式进行迭代求解得到散射估计值。
在本公开的一实施例中，所述能谱校正后探测器像素探测得到的总光子数的表达式为：能谱校正后探测器像素探测得到的总光子数为能谱校正后的透射光子数与能谱校正后散射值分布函数的叠加。
在本公开的一实施例中，对应于两个不同位置的光源的情况，所述能谱校正后探测器像素探测得到的总光子数的表达式满足：
其中，
表示当光源Pi闪烁时，经过能谱校正后的探测器像素探测得到的总光子数，i表示不同位置光源的序号，个数为至少两个，这里i＝1，2；
表示光源Pi发光时穿过滤波栅格的光子数；fi(·)为光源Pi闪烁时对应的能谱校正函数；
表示光源Pi闪烁时实际测算得到的投影值；
表示光源Pi闪烁时对应的透射投影值；Gi(Is)表示光源Pi闪烁时经过能谱校正之后关于散射值IS的分布函数，该分布函数的表达式有多种，根据能谱硬化对于散射分布的影响确定；
对应于个数大于两个的不同位置的光源的情况，按照上述总光子数的表达式类推增加方程的个数。
在本公开的一实施例中，所述基于预先建立的投影数据值与散射值映射关系查找表和测量的投影值得到散射估计值，包括：
预先建立投影数据值与散射值之间的映射关系查找表：将被扫描物体放于多光源CT成像系统中，获取至少两个不同位置的光源在不设置散射遮挡硬件时，不放置和放置被扫描物体对应在探测器同一像素点处的投影数据值，不放置被扫描物体对应的投影数据值为光源Pi发光时穿过滤波栅格的光子数
放置被扫描物体对应的投影数据值为探测器像素探测得到的总光子数
i表示不同位置光源的序号，个数为至少两个；获取在设置散射遮挡硬件时放置被扫描物体对应在探测器同一像素点的透射光子数
根据不设置散射遮挡硬件时光源发光时穿过滤波栅格的光子数
不设置散射遮挡硬件时探测器像素探测得到的总光子数
和设置散射遮挡硬件时同一探测器像素探测得到的透射光子数
得到经过能谱校正后的探测器像素探测得到的总光子数
和透射光子数
利用经过能谱校正后的探测器像素探测得到的总光子数
减去经过能谱校正后的探测器像素探测得到的透射光子数
求得散射值，该散射值为经过能谱校正后的散射值；建立探测器像素探测得到的总光子数
和所述散射值之间的映射关系查找表；以及
将被扫描物体放置于多光源CT成像系统中，实际测量至少两个不同位置的光源对应探测器同一像素点的投影值，并基于所述查找表根据测量得到的投影值得到该投影值对应的散射估计值。
在本公开的一实施例中，所述经过能谱校正后散射值的分布函数基于能谱硬化对于散射分布的影响确定，包括：在能谱硬化对于散射分布具有影响的情况下，则假设经过能谱校正后散射值的分布函数具有变化，将实际测量投影值表示为透射投影值减去散射因素引起的偏差的情况下，将能谱校正后实际测量投影值采用一阶泰勒展开，建立能谱校正后实际测量投影值与能谱校正后透射投影值间的方程，将能谱校正后实际测量总光子数的表达式与所述方程联立进而得到能谱校正后的散射分布函数。
在本公开的一实施例中，所述经过能谱校正后散射值的分布函数基于能谱硬化对于散射分布的影响确定，包括：在能谱硬化对于散射分布无影响的情况下，则假设经过能谱校正后散射值的分布函数无变化，即得到如下关系：能谱校正后散射值的分布函数为初始散射值。
在本公开的一实施例中，所述根据能谱校正后散射值的表达式进行迭代求解得到散射估计值，包括：在每个迭代循环内，基于能谱校正后散射值的表达式获取上一步得到的累积散射，对于第一步迭代，上一步得到的累积散射为散射值初值，每一步迭代利用上一步得到的累积散射计算散射去除后的投影值，并通过该投影值计算当前步残留的散射，最后更新当前步的累积散射，直到到达给定的迭代次数或者当前步计算的残留散射小于给定值时停止迭代。
在本公开的一实施例中，所述对至少两个不同位置的光源在探测器同一像素点处的投影数据进行能谱校正的步骤中，基于如下基本假设：处于不同位置的光源分别发光时，对于同一个探测器像素而言，射线经过被扫描物体的路径长度相等，同时处于不同位置的光源对应的散射光子数相等。
在本公开的一实施例中，该CT散射校正方法还包括：利用求解得到的散射估计值对投影数据进行校正，得到经过散射校正后的投影值。
在本公开的一实施例中，所述将被扫描物体放于多光源CT成像系统中，获取至少两个不同位置的光源在探测器同一像素点处的投影数据，包括：将被扫描物体放于多光源CT成像系统中，不同位置光源发射的X射线经过非均匀滤波栅格进行能谱调制，能谱调制后的X射线经过被扫描物体吸收之后由探测器进行探测，分别获取至少两个不同位置光源对应在探测器同一像素点处的投影数据。
根据本公开的另一个方面，提供了一种CT散射校正系统，用于执行本公开提及的任一种CT散射校正方法，所述CT散射校正系统包括：多光源CT成像系统，用于放置被扫描物体，以获取至少两个不同位置的光源在探测器同一像素点处的投影数据；能谱校正模块，用于对至少两个不同位置的光源在探测器同一像素点处的投影数据进行能谱校正，获得能谱校正后探测器像素探测得到的总光子数的表达式，该表达式中包括：经过能谱校正后散射值的分布函数，其中所述经过能谱校正后散射值的分布函数基于能谱硬化对于散射分布的影响确定；以及散射估计值求解模块，用于根据能谱校正后散射值的分布函数求解出能谱校正后散射值的表达式，并根据能谱校正后散射值的表达式进行迭代求解得到散射估计值；和/或，所述CT散射校正系统包括：存储模块，用于存储预先建立的投影数据值与散射值之间的映射关系查找表，映射关系中的所述散射值为经过能谱校正后的散射值；多光源CT成像系统，用于放置被扫描物体，以获取至少两个不同位置的光源对应探测器同一像素点的投影值；以及散射估计值计算模块，用于基于所述查找表根据测量得到的投影值得到该投影值对应的散射估计值。
在本公开的一实施例中，所述CT散射校正系统还包括：投影数据校正模块，用于利用求解得到的散射估计值对投影数据进行校正，得到经过散射校正后的投影值。
在本公开的一实施例中，所述多光源CT成像系统包括：多个光源、非均匀滤波栅格以及探测器，将被扫描物体放置于非均匀滤波栅格与探测器之间；其中，所述多个光源为分布式射线源，该分布式射线源中各个射线源按照设定依次闪烁；或者所述多个光源由单个射线源的焦点快速移动形成。
(三)有益效果
从上述技术方案可以看出，本公开提供的CT散射校正方法及系统，具有以下有益效果：
提出两种方式对散射值进行估计，一种方式是基于预先建立的投影数据值与散射值映射关系查找表和测量的投影值得到散射估计值，映射关系中的所述散射值为经过能谱校正后的散射值；另一种方式是基于能谱物理对散射值的分布函数进行校正后迭代求解得到散射估计值，通过获取至少两个不同位置的光源在探测器同一像素点处的投影数据，对至少两个不同位置的光源在探测器同一像素点处的投影数据进行能谱校正，校正后的表达式中涵盖了能谱硬化对于散射分布的影响，根据能谱硬化对于散射分布的影响可以确定经过能谱校正后散射值的分布函数，进而根据能谱校正后散射值的分布函数求出能谱校正后散射值的表达式，对能谱校正后散射值的表达式进行迭代求解可以得到散射估计值。两种方法均通过在散射估计过程中加入能谱校正的因素，给出了能谱硬化对于散射分布影响的计算表达形式或者测量得到方式，基于不同的能谱硬化对于散射分布影响的假设，按照迭代求解的方式或者基于散射遮挡硬件测量的方式可以得到能谱校正后的散射估计值，对于硬件的要求不高，同时还具有迭代效率高，计算准确度较高的优点。
附图说明
图1为现有技术中常用CT系统结构及散射分布示意图。
图2为根据本公开一实施例所示的多光源CT成像系统的结构示意图。
图3为根据本公开一实施例所示的CT散射校正方法流程图。
图4为根据本公开一实施例所示的CT散射校正方法还包括预先获取能谱校正函数的步骤的流程图。
图5为根据本公开一实施例所示的CT散射校正系统的结构框图。
【符号说明】
2-CT校正系统；
21-多光源CT成像系统；
22a-能谱校正模块；
22b-存储模块；
23-散射估计值求解模块；
24-投影数据校正模块。
具体实施方式
图1为现有技术中常用CT系统结构及散射分布示意图，参照图1所示，进行CT成像时，由射线源发出X射线，X射线经过被扫描物体的衰减之后被探测器探测得到投影数据。X射线在经过被扫描物体时有部分光子发生透射，部分光子发生散射，散射分布相比于通常穿过物体进入探测器的投影数据分布，变化比较缓慢，因此以高频主束分布来表示透射信息，以低频散射分布来表示散射信息，由探测器测试得到的信息中包括高频主束分布和低频散射分布。其中，发生透射的透射光子为有用信息，透射光子连同散射光子一同穿过待测物体进入探测器，散热光子会造成图像伪影和CT值不准确等问题。
目前X射线作为CT成像的射线源，由X射线管发射的X射线是分布在频谱上的量，并非严格的单值。当X射线通过吸收物质时，X射线频谱的不同部分被吸收(或称为衰减)的量并不相同，较低能量的X射线优先被吸收，对应X射线频谱分布的峰值会向高能方向移动，这一过程称之为能谱硬化。射线束硬化使得X射线对应的光子吸收不均衡，这种非线性衰减不作修正，会产生环状或条状伪影，由此导致对CT成像的质量造成很大影响。
本公开提供了一种CT散射校正方法及系统，考虑到能谱硬化和X射线散射引起的图像伪影和CT值不准确的问题，将被扫描物体放于多光源CT成像系统中，不同位置光源发射的X射线经过非均匀滤波栅格进行能谱调制，能谱调制后的X射线经过被扫描物体吸收之后由探测器进行探测，一种途径为采用查找表的方法，预先建立投影数据值与散射值之间的映射关系查找表；以及将被扫描物体放置于多光源CT成像系统中，实际测量至少两个不同位置的光源对应探测器同一像素点的投影值，并基于所述查找表根据测量得到的投影值得到该投影值对应的散射估计值。预先建立投影数据值与散射值之间的映射关系查找表的方法包括：将被扫描物体放于多光源CT成像系统中，获取至少两个不同位置的光源在不设置散射遮挡硬件时，不放置和放置被扫描物体对应在探测器同一像素点处的投影数据值，不放置被扫描物体对应的投影数据值为光源发光时穿过滤波栅格的光子数，放置被扫描物体对应的投影数据值为探测器像素探测得到的总光子数；获取在设置散射遮挡硬件时放置被扫描物体对应在探测器同一像素点的透射光子数，根据不设置散射遮挡硬件时光源发光时穿过滤波栅格的光子数、不设置散射遮挡硬件时探测器像素探测得到的总光子数和设置散射遮挡硬件时同一探测器像素探测得到的透射光子数得到经过能谱校正后的探测器像素探测得到的总光子数和经过能谱校正后的探测器像素探测到的透射光子数，利用经过能谱校正后的探测器像素探测得到的总光子数减去经过能谱校正后的探测器像素探测得到的透射光子数求得散射值，该散射值为经过能谱校正后的散射值；建立探测器像素探测得到的总光子数和所述散射值之间的映射关系查找表。
另外一种途径为采用基于能谱物理的散射估计方法，将被扫描物体放于多光源CT成像系统中，获取至少两个不同位置的光源在探测器同一像素点处的投影数据；对至少两个不同位置的光源在探测器同一像素点处的投影数据进行能谱校正，获得能谱校正后探测器像素探测得到的总光子数的表达式，该表达式中包括：经过能谱校正后散射值的分布函数，其中所述经过能谱校正后散射值的分布函数基于能谱硬化对于散射分布的影响确定；以及根据能谱校正后散射值的分布函数求解出能谱校正后散射值的表达式，并根据能谱校正后散射值的表达式进行迭代求解得到散射估计值。上述CT散射校正方法中，不局限于两个不同位置光源对应的情况，还可以拓展至N个不同位置光源的情况进行散射校正。此种散射估计方法对应的迭代次数较少，迭代效率很高。当然，在下面的实施例中，对于能谱校正后散射值的分布函数仅提供了两种示例性的实例，本公开的能谱校正后散射值的分布函数Gi(Is)的表达式不局限于列举的实施例，在本公开的其他实施方式中，可以对能谱校正后散射值的分布函数表示为其他表达式，进行散射值估计的过程与示例的方式一致，仅需要更换能谱校正后散射值的分布函数Gi(Is)以及散射值IS的表达式即可，这里不再穷举。
为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开进一步详细说明。
第一实施例
在本公开的第一个示例性实施例中，提供了一种CT散射校正方法。
图2为根据本公开一实施例所示的多光源CT成像系统的结构示意图。图3为根据本公开一实施例所示的CT散射校正方法流程图。
本实施例中，基于多光源CT成像系统来实施上述CT散射校正方法。参照图2所示，该多光源CT成像系统包括：多个光源、非均匀滤波栅格以及探测器。在CT成像过程中，X射线源与探测器围绕物体进行相对旋转运动，以此获得不同旋转角度下的CT数据，即为CT投影值。不同位置光源发射的X射线经过非均匀滤波栅格进行能谱调制，能谱调制后的X射线经过被扫描物体吸收之后由探测器进行探测。
本实施例中，将被扫描物体放于多光源CT成像系统中，例如图2中所示意的将被扫描物体放置于非均匀滤波栅格与平板探测器之间。多光源扫描的一种实现方式为：多个光源为分布式射线源，该分布式射线源中各个射线源依次闪烁，射线源可为X光管、碳纳米管或加速器，光源分布形式及分布间距可以根据需求进行定义。多光源扫描的另一种实现方式为：单个射线源的焦点快速移动，以此改变射线源发射射线的初始位置，焦点可以根据需求自定义移动方向。
非均匀滤波栅格为非均匀排布的滤波栅格，滤波栅格由能够吸收部分射线的材料制造，能够对入射的射线能谱进行调制，根据与射线源焦点的相对位置产生不同能谱的射线，实现对于入射射线的能谱调制。滤波栅格的材料包括但不限于为如下材料的一种或几种：铜、银、锡以及钼等。
参照图3所示，本实施例的CT散射校正方法，包括：步骤S111、步骤S112、步骤S113和步骤S114。
在步骤S111，将被扫描物体放于多光源CT成像系统中，获取至少两个不同位置的光源在探测器同一像素点处的投影数据。下面以两个不同位置的光源对应的情况作为示例，需要说明的是，不同位置的光源的个数可以拓展至多个(大于2个)。
参照图2所示，将被扫描物体放于多光源CT成像系统中，不同位置光源发射的X射线经过非均匀滤波栅格进行能谱调制，能谱调制后的X射线经过被扫描物体吸收之后由探测器进行探测，分别获取至少两个不同位置光源对应在探测器同一像素点处的投影数据。
本实施例中，在同一个探测器像素点可以获取至少两个不同位置的光源下的X光信号，该探测器像素获取的X光信号即为投影数据。该投影数据中包含透射光子的贡献和散射光子的贡献，在后续的步骤中需要去除散射光子的部分，以校正投影数据。
参照图2所示，该多光源CT成像系统中，示意了相邻的光源Pn和光源Pn+1经过非均匀滤波栅格和被扫描物体后入射到探测器的路径，n为正整数，其他任意两个不同位置的光源Pa和Pb的情况与上述光源Pn和光源Pn+1的情况类似，这里仅以光源Pn和光源Pn+1这两个光源作为示例。当两个任意位置的光源分别闪烁时，例如光源P1和光源P2分别闪烁，对于同一个探测器像素，采集得到的X光信号包括透射光子的贡献以及散射光子的贡献，可以简单地表示为如下表达式：
其中，
表示当光源Pi闪烁时，探测器像素探测得到的总光子数，其中i＝1，2；
表示当光源Pi闪烁时对应的散射光子数；
表示当光源Pi闪烁时对应的透射光子数。
透射光子数可以进一步表示为如下表达式：
其中，Si(E)表示光源Pi照向探测器的入射能谱，um(E)表示滤波栅格在能量E下的线性衰减系数；
表示射线路径上滤波栅格的厚度；uo(E)表示被扫描物体在能量E下的线性衰减系数；Li表示射线穿过被扫描物体的路径长度。考虑到CT成像系统中不同方向上的出射能谱差异性以及不同探测器之间的不一致性，作为优选，公式(2)中的入射能谱可以是等效入射能谱，通过对每个探测器像素探测得到的初始能谱进行优化而得到该等效入射能谱。另外，光源在不同位置下射线所经过滤波栅格的厚度可以通过事先标定和优化得到。
上述公式(2)中的透射光子数的表达式中包含滤波栅格和被扫描物体的信息，后续步骤旨在将散射光子数进行能谱校正，从而对应得到去除了散射光子贡献的投影数据，即得到相对准确的透射光子数，从而可以根据公式(2)来进行CT成像的过程。
术语“入射能谱”是指相对于探测器为入射的方向，由探测器探测到的入射能谱；相对于光源来说，从光源发出。
在步骤S112，对至少两个不同位置的光源在探测器同一像素点处的投影数据进行能谱校正，获得能谱校正后探测器像素探测得到的总光子数的表达式，该表达式中包括：经过能谱校正后散射值的分布函数。
本实施例中，进行能谱校正后的表达式将探测器像素探测得到的总光子数近似表示为两部分的加和，一部分为光源发光时穿过滤波栅格的光子数与能谱校正因子的乘积；另一部分为经过能谱校正后散射值的分布函数。
对于同一个探测器像素而言，当来自不同位置光源的经过扫描物体的射线路径在空间上比较接近时，散射基本相同：
L1＝L2＝L， (3)
IS1＝IS2＝IS， (4)
其中，L为射线穿过被扫描物体的路径长度值；IS为散射值。
以上可以通过合理配置不同光源的间距实现。此时可以认为当光源P1和光源P2分别发光时，对于同一个探测器像素而言，射线经过被扫描物体的路径长度基本相等，不同光源闪烁时对应的散射光子数也相等。
将公式(4)带入至公式(1)中，得到如下表达式：
考虑到能谱硬化对于散射分布的影响，利用能谱校正的方法对探测得到的光子数进行能谱校正。能谱校正后探测器像素探测得到的总光子数的表达式为：能谱校正后探测器像素探测得到的总光子数为能谱校正后的透射光子数与能谱校正后散射值分布函数的叠加。
在一实施例中，对公式(5)进行能谱校正之后，得到如下表达式：
其中，
表示当光源Pi闪烁时，经过能谱校正后的探测器像素探测得到的总光子数，i＝1，2；
表示光源Pi发光时穿过滤波栅格的光子数；fi(·)为光源Pi闪烁时对应的能谱校正函数，可以由预先模拟或实验标定得到；
表示光源Pi闪烁时实际测算得到的投影值；
表示光源Pi闪烁时对应的透射投影值；Gi(Is)表示光源Pi闪烁时经过能谱校正之后关于散射值IS的分布函数。上述约等号可以视为等号，表示进行能谱校正的过程。
图4为根据本公开一实施例所示的CT散射校正方法还包括预先获取能谱校正函数的步骤的流程图。
需要说明的是，上述能谱校正函数fi(·)可以由预先模拟或实验标定得到，例如参照图4所示，可以在步骤S101中，基于入射能谱(优选等效入射能谱)和滤波栅格厚度预先得到能谱校正函数。
其中，
和
分别满足如下表达式：
在步骤S113，对经过能谱校正后散射值的分布函数的变化情况进行假设，根据对应的假设推导得出能谱校正后散射值的表达式，并根据能谱校正后散射值的表达式进行迭代算法求解得出散射估计值。该步骤S113中，由于散射值的分布函数已经经过了能谱校正，则根据对应假设推导得出的散射值已经是经过能谱校正后的散射值，由该散射值最终迭代求解得到的散射估计值已经涵盖了能谱硬化的影响。
基于Gi(Is)取值的不同，可以进行不同的迭代求解过程。
下面采用两个实施例来示例性说明两种迭代方法对散射值进行估计的过程。其他形式的Gi(Is)的表达式对应的求解过程与下面的示例类似，只需要代入相应Gi(Is)和Is的表达式中即可。该步骤S113中，对经过能谱校正后散射值的分布函数的变化情况进行假设，一种假设是能谱硬化对于散射分布没有影响或者影响很小，即经过能谱校正后散射值的分布函数不发生变化。另一种假设是能谱硬化对于散射分布具有影响，即经过能谱校正后散射值的分布函数发生变化。
下面针对不同的假设来详细介绍如何根据对应的假设推导得出散射值的表达式，并根据散射值的表达式对散射值进行迭代算法求解得出散射值估计的过程。根据前述能谱硬化对于散射分布的影响得出的能谱校正后散射值的分布函数Gi(Is)，可以求解得出由投影数据计算散射的公式，即得到散射值的表达式，该散射值为经过能谱校正后的散射值，然后基于散射值的表达式可以通过迭代算法求解得出散射估计值。
第一种迭代方法示例如下，所述经过能谱校正后散射值的分布函数基于能谱硬化对于散射分布的影响确定，包括：在能谱硬化对于散射分布无影响的情况下，则假设经过能谱校正后散射值的分布函数无变化，即得到如下关系：能谱校正后散射值的分布函数不发生变化，该能谱校正后散射值的分布函数表示为初始散射值，这里的初始散射值含义为并未发生变化的散射值Is，在公式(10)中以Is表示初始散射值。
该实施例中，可以简单地认为经过能谱校正之后散射的分布没有改变，即存在如下关系式：
Gi(Is)＝Is， (10)
则将公式(10)代入至公式(6)中，得到如下表达式：
对于同一个探测器像素而言，两个不同位置的光源出射的射线经过扫描物体的路径长度基本一致，那么经过能谱校正函数会映射到相同的投影值，即存在如下关系式：
那么根据公式(6)、(11)和(12)可以求解得到能谱校正后Is的表达式如下：
针对公式(13)所示的能谱校正后Is的表达式，采用迭代算法进行求解。
在本公开的一实施例中，根据能谱校正后散射值的表达式进行迭代求解得到散射估计值，包括：在每个迭代循环内，基于能谱校正后散射值的表达式获取上一步得到的累积散射，对于第一步迭代，上一步得到的累积散射为散射值初值，每一步迭代利用上一步得到的累积散射计算散射去除后的投影值，并通过该投影值计算当前步残留的散射，最后更新当前步的累积散射，直到到达给定的迭代次数或者当前步计算的残留散射小于给定值时停止迭代。
下面以具体迭代的实例进行说明。在一实施例中，根据散射值的表达式对散射值进行迭代算法求解得出散射估计值，包括：
给散射值赋予初值Is 0，初始化设定停止准则stop_c以及最大迭代步数Niter；迭代步数从1开始，每次循环增加1步，从1增加至Niter，即满足：迭代步数iter＝1：Niter，这里的等号表示赋值；
在每个迭代循环内，将特定迭代步数下的散射值代入至
的表达式中得到该特定迭代步数下
的值，将
的值代入至所述散射值的表达式中得到特定迭代步数下的散射值残差ΔIs iter；则根据得到的特定迭代步数下的散射值残差ΔIs iter与上一步迭代得到的散射值
计算得到本次迭代循环的散射值
本次迭代循环的散射值
满足：
如果在某步迭代后，|ΔIs|小于停止准则stop_c，则停止迭代，得到最终的散射估计值为
否则，一直迭代至最大迭代步数，返回最终的散射估计值为
下面采用伪代码的形式进行示例性说明。在一实例中，迭代算法求解散射值Is的伪代码(Pseudocode)如下：
(a1)初始化Is 0、停止准则stop_c以及最大迭代步数Niter；
(b1)for iter＝1：Niter：
(c1)
(d1)
(e1)
(f1)
(g1)如果|ΔIs|小于stop-c，则停止迭代，最后的散射估计值为
(h1)endfor.返回最后的散射估计值为
第二种迭代方法示例如下，在能谱硬化对于散射分布具有影响的情况下，则假设经过能谱校正后散射值的分布函数具有变化，将实际测量投影值表示为透射投影值减去散射因素引起的偏差的情况下，将能谱校正后实际测量投影值采用一阶泰勒展开，建立能谱校正后实际测量投影值与能谱校正后透射投影值间的方程，将能谱校正后实际测量总光子数的表达式与所述方程联立进而得到能谱校正后的散射分布函数。
具体而言，该实施例中，考虑到能谱硬化对于散射分布的影响，并对散射进行相应的补偿。透射投影值
和实际测量投影值
之间存在散射因素引起的偏差Δp，利用公式表示如下：
对于实际测量投影值
进行能谱校正，并应用一阶泰勒展开进行近似，能谱校正后实际测量投影值采用一阶泰勒展开的表达式如下：
其中，f′表示f的一阶导数。
在上述(16)表达式的基础上，可以建立能谱校正后实际测量投影值与能谱校正后透射投影值间的方程。具体建立过程为：根据公式(16)可以进行推导，得到能谱校正之后的探测器像素探测得到的总光子数
的表达式如下：
基于公式(17)中的泰勒近似，对应公式(6)可以得到如下表达式：
方程(17)示例性描述了本实施例中能谱校正后实际测量投影值与能谱校正后透射投影值间的方程，其中本领域技术人员公知，实际测量投影值与实际测量时探测器像素探测得到的总光子数之间存在负对数关系。将能谱校正后实际测量总光子数对应的表达式(6)与所述方程(17)联立进而得到能谱校正后散射值的分布函数(散射值的分布函数也可以称为散射分布函数)，即得到公式(18)。
然后基于能谱校正后散射值的分布函数求解出散射值的表达式。在本实施例中，可以将不同像素点对应的能谱校正后散射值的分布函数代入能谱校正后实际测量时探测器像素探测得到的总光子数的表达式中，从而求解得出能谱校正后散射值的表达式。具体而言，将公式(18)代入公式(6)中，可以得到如下表达式：
根据公式(19)列出的方程可以解出Is的表达式。下面列出求解出的能谱校正后Is的表达式如下：
针对公式(20)所示的能谱校正后Is的表达式，采用迭代算法进行求解。
根据能谱校正后散射值的表达式进行迭代求解得到散射估计值，包括：在每个迭代循环内，基于能谱校正后散射值的表达式获取上一步得到的累积散射，对于第一步迭代，上一步得到的累积散射为散射值初值，每一步迭代利用上一步得到的累积散射计算散射去除后的投影值，并通过该投影值计算当前步残留的散射，最后更新当前步的累积散射，直到到达给定的迭代次数或者当前步计算的残留散射小于给定值时停止迭代。
下面以具体实例进行介绍。在一实例中，迭代算法求解散射值Is的伪代码(Pseudocode)如下：
(a2)初始化Is 0、停止准则stop_c以及最大迭代步数Niter；
(h2)for iter＝1：Niter.
(c2)
(d2)
(e2)
(f2)
(g2)如果|ΔIs|小于stop_c，则停止迭代，最后的散射估计值为
(h2)end for.返回最后的散射估计值为
上述实施例中，第二种迭代方法具有更高的迭代效率，可以很快逼近散射估计值。
在步骤S114，利用求解得到的散射估计值对投影数据进行校正，得到经过散射校正后的投影值。
本实施例中，步骤S114中的投影数据中包含透射光子的贡献和散射光子的贡献，经过散射校正后的投影值仅包含透射光子的贡献。其中利用求解得到的散射估计值对投影数据校正，包括：在投影数据中减去散射估计值，以实现对于投影数据的校正。
当然，上述实施例中，以两个不同位置的光源在探测器同一像素点处的投影数据列出两个方程进行示例，本公开中，i表示不同位置光源的序号，个数为至少两个，这里仅以i＝1，2作为示例，实际上，对于个数大于两个的情况，i＝1，2，…，k，可以按照公式(6)中的形式进行类推得到包含k个方程的方程组，k＞2；
在不同位置的光源的个数大于2的情况下，不仅可以按照上述方法求解得到散射估计值，还可以得出其他参数。
第二实施例
本公开的第二个示例性实施例中，提供了一种CT散射校正方法，本实施例的CT散射校正方法与第一实施例不同，是利用查找表法来对散射进行校正。
本实施例中，基于预先建立的投影数据值与散射值映射关系查找表和测量的投影值得到散射估计值，包括：预先建立投影数据值与散射值之间的映射关系查找表；以及实际测量被扫描物体放置于多光源CT成像系统中，至少两个不同位置的光源对应探测器同一像素点的投影值，并基于所述查找表根据测量得到的投影值得到该投影值对应的散射估计值。
在一实例中，预先建立投影数据值与散射值之间的映射关系查找表，包括：将被扫描物体放于多光源CT成像系统中，获取至少两个不同位置的光源在不设置散射遮挡硬件时，不放置和放置被扫描物体对应在探测器同一像素点处的投影数据值。其中不放置被扫描物体对应的投影数据值为光源Pi发光时穿过滤波栅格的光子数
放置被扫描物体对应的投影数据值为探测器像素探测得到的总光子数
i表示不同位置光源的序号，个数为至少两个；获取在设置散射遮挡硬件时放置被扫描物体对应在探测器同一像素点的透射光子数
根据
和
得到经过能谱校正后的探测器像素探测得到的总光子数
和经过能谱校正后的探测器像素探测得到的透射光子数
可以参照前述公式(6)的求得。利用经过能谱校正后的探测器像素探测得到的总光子数
减去经过能谱校正后的探测器像素探测得到的透射光子数
求得散射值，该散射值为经过能谱校正后的散射值。基于上述，可以建立探测器像素探测得到的总光子数
和能谱校正后散射值之间的映射关系查找表。也就是说，建立了不设置散射遮挡硬件、放置被扫描物体时对应的实测投影数据值和能谱校正后散射值之间的查找表，直接可以利用实际测量的投影值和查找表获取散射估计值，该散射估计值为能谱校正后的散射估计值，无需重复进行测量或者其他复杂计算过程。
然后，将被扫描物体放置于多光源CT成像系统中，实际测量至少两个不同位置的光源对应探测器同一像素点的投影值，并基于所述查找表根据测量得到的投影值得到该投影值对应的散射估计值。
第三实施例
本公开的第三个示例性实施例中，提供了一种CT散射校正系统，用于执行本公开提及的任一种CT散射校正方法。
图5为根据本公开一实施例所示的CT散射校正系统的结构框图。
参照图5所示，本实施例的CT散射校正系统2包括：多光源CT成像系统21、能谱校正模块22a和散射估计值求解模块23。
其中，多光源CT成像系统21用于放置被扫描物体，以获取至少两个不同位置的光源在探测器同一像素点处的投影数据。能谱校正模块22a用于对至少两个不同位置的光源在探测器同一像素点处的投影数据进行能谱校正，获得能谱校正后探测器像素探测得到的总光子数的表达式。该表达式中包括：经过能谱校正后散射值的分布函数。其中，经过能谱校正后散射值的分布函数基于能谱硬化对于散射分布的影响确定。散射估计值求解模块23用于根据能谱校正后散射值的分布函数求解出能谱校正后散射值的表达式，并根据能谱校正后散射值的表达式进行迭代求解得到散射估计值。
在另一实施例中，CT散射校正系统2可以用于执行查找表法来实现散射值的估计校正。该CT散射校正系统2还可以包括：存储模块22b，该存储模块22b用于存储预先建立的投影数据值与散射值之间的映射关系查找表，映射关系中的所述散射值为经过能谱校正后的散射值。上述散射估计值计算模块23还用于基于所述查找表根据测量得到的投影值得到该投影值对应的散射估计值。当然，该CT散射校正系统可以同时包括：多光源CT成像系统21、能谱校正模块22a、散射估计值求解模块23和存储模块22b，或者，可以包括：多光源CT成像系统21、能谱校爪模块22a和散射估计值求解模块23或者包括：多光源CT成像系统21、存储模块22b和散射估计值求解模块23。
在本公开的一实施例中，参照图5所示，CT散射校正系统2还包括：投影数据校正模块24。该投影数据校正模块24用于利用求解得到的散射估计值对投影数据进行校正，得到经过散射校正后的投影值。
在本公开的一实施例中，参照图1所示，多光源CT成像系统21包括：多个光源、非均匀滤波栅格以及探测器，将被扫描物体放置于非均匀滤波栅格与探测器之间。本实施例中，多个光源可以为分布式射线源，该分布式射线源中各个射线源按照设定依次闪烁；或者所述多个光源可以由单个射线源的焦点快速移动形成。
综上所述，本公开提供了一种CT散射校正方法及系统，提供了两种途径来进行散射分布的估计，一种是通过查找表的方法，另一种是通过基于能谱物理的散射估计方法，在基于能谱物理的散射估计方法中，通过获取至少两个不同位置的光源在探测器同一像素点处的投影数据，对至少两个不同位置的光源在探测器同一像素点处的投影数据进行能谱校正，校正后的表达式中涵盖了能谱硬化对于散射的影响，根据能谱硬化对于散射分布的影响可以确定经过能谱校正后散射值的分布函数，进而根据能谱校正后散射值的分布函数求解出能谱校正后散射值的表达式，并根据能谱校正后散射值的表达式进行迭代求解可以得到散射估计值，两种方法均通过在散射估计过程中加入能谱校正的因素，给出了能谱硬化对于散射分布影响的计算表达形式或者测量得到方式，基于不同的能谱硬化对于散射分布影响的假设，按照迭代求解的方式或者基于散射遮挡硬件测量的方式可以得到能谱校正后的散射估计值，对于硬件的要求不高，同时还具有迭代效率高，计算准确度较高的优点。
本领域那些技术人员可以理解，可以对实施例中的系统中的模块进行自适应性地改变并且把它们设置在与该实施例不同的一个或多个模块中。可以把实施例中的模块组合成一个模块，以及此外可以把它们分成多个子模块。除了这样的特征和/或过程或者单元中的至少一些是相互排斥之外，可以采用任何组合对本说明书(包括伴随的权利要求、摘要和附图)中公开的所有特征以及如此公开的任何方法或者装置的所有过程或单元进行组合。除非另外明确陈述，本说明书(包括伴随的权利要求、摘要和附图)中公开的每个特征可以由提供相同、等同或相似目的的替代特征来代替。根据本公开实施例的模块中的任意一个或多个可以至少被部分地实现为硬件电路，例如现场可编程门阵列(FPGA)、可编程逻辑阵列(PLA)、片上系统、基板上的系统、封装上的系统、专用集成电路(ASIC)，或可以通过对电路进行集成或封装的任何其他的合理方式的硬件或固件来实现，或以软件、硬件以及固件三种实现方式中任意一种或以其中任意几种的适当组合来实现。或者，根据本公开实施例的模块中的一个或多个可以至少被部分地实现为计算机程序模块，当该计算机程序模块被运行时，可以执行相应的功能。
再者，单词“包含”或“包括”不排除存在未列在权利要求中的元件或步骤。位于元件之前的单词“一”或“一个”不排除存在多个这样的元件。
除非存在技术障碍或矛盾，本发明的上述各种实施方式可以自由组合以形成另外的实施例，这些另外的实施例均在本发明的保护范围中。
以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本公开的具体实施例而已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。
Claims (9)
1.一种CT散射校正方法，其特征在于，包括：
基于能谱物理对散射值的分布函数进行校正后迭代求解得到散射估计值；
其中，所述基于能谱物理对散射值的分布函数进行校正后迭代求解得到散射估计值，包括：
将被扫描物体放于多光源CT成像系统中，获取至少两个不同位置的光源在探测器同一像素点处的投影数据；
对至少两个不同位置的光源在探测器同一像素点处的投影数据进行能谱校正，获得能谱校正后探测器像素探测得到的总光子数的表达式，该表达式中包括：经过能谱校正后散射值的分布函数，其中所述经过能谱校正后散射值的分布函数基于能谱硬化对于散射分布的影响确定；以及
根据能谱校正后散射值的分布函数求解出能谱校正后散射值的表达式，并根据能谱校正后散射值的表达式进行迭代求解得到散射估计值。
2.根据权利要求1所述的CT散射校正方法，其特征在于，所述能谱校正后探测器像素探测得到的总光子数的表达式为：能谱校正后探测器像素探测得到的总光子数为能谱校正后的透射光子数与能谱校正后散射值分布函数的叠加。
3.根据权利要求1或2所述的CT散射校正方法，其特征在于，所述经过能谱校正后散射值的分布函数基于能谱硬化对于散射分布的影响确定，包括：在能谱硬化对于散射分布具有影响的情况下，则假设经过能谱校正后散射值的分布函数具有变化，将实际测量投影值表示为透射投影值减去散射因素引起的偏差的情况下，将能谱校正后实际测量投影值采用一阶泰勒展开，建立能谱校正后实际测量投影值与能谱校正后透射投影值间的方程，将能谱校正后实际测量总光子数的表达式与所述方程联立进而得到能谱校正后的散射分布函数。
4.根据权利要求1或2所述的CT散射校正方法，其特征在于，所述经过能谱校正后散射值的分布函数基于能谱硬化对于散射分布的影响确定，包括：在能谱硬化对于散射分布无影响的情况下，则假设经过能谱校正后散射值的分布函数无变化，即得到如下关系：能谱校正后散射值的分布函数为初始散射值。
5.根据权利要求1或2所述的CT散射校正方法，其特征在于，所述根据能谱校正后散射值的表达式进行迭代求解得到散射估计值，包括：在每个迭代循环内，基于能谱校正后散射值的表达式获取上一步得到的累积散射，对于第一步迭代，上一步得到的累积散射为散射值初值，每一步迭代利用上一步得到的累积散射计算散射去除后的投影值，并通过该投影值计算当前步残留的散射，最后更新当前步的累积散射，直到到达给定的迭代次数或者当前步计算的残留散射小于给定值时停止迭代。
6.根据权利要求1所述的CT散射校正方法，其特征在于，还包括：利用求解得到的散射估计值对投影数据进行校正，得到经过散射校正后的投影值。
7.根据权利要求1所述的CT散射校正方法，其特征在于，所述将被扫描物体放于多光源CT成像系统中，获取至少两个不同位置的光源在探测器同一像素点处的投影数据，包括：
将被扫描物体放于多光源CT成像系统中，不同位置光源发射的X射线经过非均匀滤波栅格进行能谱调制，能谱调制后的X射线经过被扫描物体吸收之后由探测器进行探测，分别获取至少两个不同位置光源对应在探测器同一像素点处的投影数据。
8.一种CT散射校正系统，其特征在于，用于执行权利要求1-7 中任一项所述的CT散射校正方法，所述CT散射校正系统包括：
多光源CT成像系统，用于放置被扫描物体，以获取至少两个不同位置的光源在探测器同一像素点处的投影数据；
能谱校正模块，用于对至少两个不同位置的光源在探测器同一像素点处的投影数据进行能谱校正，获得能谱校正后探测器像素探测得到的总光子数的表达式，该表达式中包括：经过能谱校正后散射值的分布函数，其中所述经过能谱校正后散射值的分布函数基于能谱硬化对于散射分布的影响确定；以及
散射估计值求解模块，用于根据能谱校正后散射值的分布函数求解出能谱校正后散射值的表达式，并根据能谱校正后散射值的表达式进行迭代求解得到散射估计值。
9.根据权利要求8所述的CT散射校正系统，其特征在于，所述多光源CT成像系统包括：多个光源、非均匀滤波栅格以及探测器，将被扫描物体放置于非均匀滤波栅格与探测器之间；其中，所述多个光源为分布式射线源，该分布式射线源中各个射线源按照设定依次闪烁；或者所述多个光源由单个射线源的焦点快速移动形成。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010133896.6A CN111643104B (zh) | 2020-02-28 | 2020-02-28 | Ct散射校正方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010133896.6A CN111643104B (zh) | 2020-02-28 | 2020-02-28 | Ct散射校正方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111643104A CN111643104A (zh) | 2020-09-11 |
| CN111643104B true CN111643104B (zh) | 2022-03-18 |
Family
ID=72346270
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010133896.6A Active CN111643104B (zh) | 2020-02-28 | 2020-02-28 | Ct散射校正方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111643104B (zh) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112415569B (zh) * | 2020-10-30 | 2022-06-17 | 湖北锐世数字医学影像科技有限公司 | 能量校正方法、装置、设备、系统及计算机可读存储介质 |
| CN113409414B (zh) * | 2021-06-08 | 2024-03-26 | 江苏一影医疗设备有限公司 | X线图像的散射修正方法、装置、电子设备、存储介质 |
| CN113237903B (zh) * | 2021-06-15 | 2022-06-10 | 清华大学 | 基于双层平板探测器锥束ct的散射校正方法及装置 |
| CN115270075B (zh) * | 2022-08-03 | 2023-04-25 | 北京朗视仪器股份有限公司 | 基于钨靶x射线球管的双能ct系统的能谱数据校正方法 |
| CN116413290B (zh) * | 2023-02-22 | 2023-11-28 | 奥影检测科技(上海)有限公司 | 一种工业ct散射校正方法 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| DE602005004653T2 (de) * | 2004-03-17 | 2009-01-29 | Philips Intellectual Property & Standards Gmbh | Korrektur der Strahlungsaufhärtung und der Dämpfung in der coherent scatter computed tomography (CSCT) |
| CN101226642A (zh) * | 2008-01-25 | 2008-07-23 | 西安交通大学 | 基于ct数据一致性的投影射束硬化校正方法 |
| US10716527B2 (en) * | 2017-02-17 | 2020-07-21 | Canon Medical Systems Corporation | Combined sinogram- and image-domain material decomposition for spectral computed tomography (CT) |
| CN109363703B (zh) * | 2018-10-18 | 2020-05-19 | 清华大学 | Ct系统能谱不一致性的校正方法 |
-
2020
- 2020-02-28 CN CN202010133896.6A patent/CN111643104B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN111643104A (zh) | 2020-09-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111643104B (zh) | Ct散射校正方法及系统 | |
| GB2551029B (en) | Method for performing material decoposition using a dual-energy X-Ray CT and corresponding Dual-energy X-Ray CT Apparatus | |
| US11331064B2 (en) | Method and apparatus for computed tomography (CT) and material decomposition with pile-up correction calibrated using real pulse pileup effect and detector response | |
| JP7139119B2 (ja) | 医用画像生成装置及び医用画像生成方法 | |
| US8705827B2 (en) | Scatter correction methods | |
| US9476993B2 (en) | Apparatus and method for computing detector response of a photon-counting detector | |
| US8958524B2 (en) | Correction of projection data in radiation system | |
| JP6386997B2 (ja) | Ｘ線コンピュータ断層撮影装置及び補正方法 | |
| US20160287205A1 (en) | Pre-reconstruction calibration, data correction, and material decomposition method and apparatus for photon-counting spectrally-resolving x-ray detectors and x-ray imaging | |
| US7496171B2 (en) | Method for estimating the radiation scattered in a two-dimensional detector | |
| US9801595B2 (en) | Count-weighted least squares parameter estimation for a photon-counting detector | |
| US10507005B2 (en) | Spectral calibration of spectral computed tomography (CT) | |
| US20200193654A1 (en) | System and method for statistical iterative reconstruction and material decomposition | |
| US20160217594A1 (en) | Cost-function based method and apparatus for projection-domain basis decomposition in spectral computed tomography | |
| US10605933B2 (en) | X-ray spectral calibration technique for cone-beam CT | |
| JP2009082616A (ja) | 多管球ｘ線ｃｔにおける散乱線強度分布のスケーリング方法および多管球ｘ線ｃｔ装置 | |
| US9230348B2 (en) | Imaging system for imaging a region of interest | |
| Dong et al. | Relationship between x‐ray illumination field size and flat field intensity and its impacts on x‐ray imaging | |
| US10646186B2 (en) | X-ray CT apparatus, information processing device and information processing method | |
| JP2016104125A (ja) | Ｘ線ｃｔ装置、画像処理装置およびプログラム | |
| US10497153B2 (en) | Heel effect correction in computed tomography | |
| US11013487B2 (en) | Method and apparatus for computed tomography (CT) and material decomposition with count-rate dependent pileup correction | |
| CN108065950B (zh) | 一种放射成像方法及其系统 | |
| US11879856B2 (en) | Two-step material decomposition calibration method for a full size photon counting computed tomography system | |
| CN111134709B (zh) | 一种多能量ct基材料分解方法 |
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 |