CN115270075B - 基于钨靶x射线球管的双能ct系统的能谱数据校正方法 - Google Patents
基于钨靶x射线球管的双能ct系统的能谱数据校正方法 Download PDFInfo
- Publication number
- CN115270075B CN115270075B CN202210927488.7A CN202210927488A CN115270075B CN 115270075 B CN115270075 B CN 115270075B CN 202210927488 A CN202210927488 A CN 202210927488A CN 115270075 B CN115270075 B CN 115270075B
- Authority
- CN
- China
- Prior art keywords
- energy spectrum
- energy
- initial
- data
- correction
- 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
- 238000001228 spectrum Methods 0.000 title claims abstract description 256
- 238000012937 correction Methods 0.000 title claims abstract description 99
- 238000000034 method Methods 0.000 title claims abstract description 62
- WFKWXMTUELFFGS-UHFFFAOYSA-N tungsten Chemical compound [W] WFKWXMTUELFFGS-UHFFFAOYSA-N 0.000 title claims abstract description 52
- 229910052721 tungsten Inorganic materials 0.000 title claims abstract description 52
- 239000010937 tungsten Substances 0.000 title claims abstract description 52
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 23
- 239000000463 material Substances 0.000 claims description 29
- 230000006870 function Effects 0.000 claims description 17
- 238000004364 calculation method Methods 0.000 claims description 12
- 239000002245 particle Substances 0.000 claims description 12
- 238000004590 computer program Methods 0.000 claims description 11
- 230000009977 dual effect Effects 0.000 claims description 10
- 238000004088 simulation Methods 0.000 claims description 7
- 230000005855 radiation Effects 0.000 description 11
- 230000003595 spectral effect Effects 0.000 description 9
- 238000000342 Monte Carlo simulation Methods 0.000 description 7
- 238000010586 diagram Methods 0.000 description 7
- 238000001914 filtration Methods 0.000 description 7
- 238000005259 measurement Methods 0.000 description 7
- 230000008569 process Effects 0.000 description 7
- 239000000126 substance Substances 0.000 description 7
- 238000004891 communication Methods 0.000 description 6
- 238000013170 computed tomography imaging Methods 0.000 description 6
- RYGMFSIKBFXOCR-UHFFFAOYSA-N Copper Chemical compound [Cu] RYGMFSIKBFXOCR-UHFFFAOYSA-N 0.000 description 5
- 238000000429 assembly Methods 0.000 description 5
- 230000000712 assembly Effects 0.000 description 5
- 229910052802 copper Inorganic materials 0.000 description 5
- 239000010949 copper Substances 0.000 description 5
- 238000012545 processing Methods 0.000 description 5
- 229910052782 aluminium Inorganic materials 0.000 description 4
- XAGFODPZIPBFFR-UHFFFAOYSA-N aluminium Chemical compound [Al] XAGFODPZIPBFFR-UHFFFAOYSA-N 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 230000005461 Bremsstrahlung Effects 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000000354 decomposition reaction Methods 0.000 description 2
- 238000003754 machining Methods 0.000 description 2
- 238000000691 measurement method Methods 0.000 description 2
- 239000007769 metal material Substances 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 239000013077 target material Substances 0.000 description 2
- 239000013076 target substance Substances 0.000 description 2
- OYPRJOBELJOOCE-UHFFFAOYSA-N Calcium Chemical compound [Ca] OYPRJOBELJOOCE-UHFFFAOYSA-N 0.000 description 1
- VVQNEPGJFQJSBK-UHFFFAOYSA-N Methyl methacrylate Chemical compound COC(=O)C(C)=C VVQNEPGJFQJSBK-UHFFFAOYSA-N 0.000 description 1
- 229920005372 Plexiglas® Polymers 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 229910052791 calcium Inorganic materials 0.000 description 1
- 239000011575 calcium Substances 0.000 description 1
- 230000001351 cycling effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 239000011521 glass Substances 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 229910052751 metal Inorganic materials 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 230000005658 nuclear physics Effects 0.000 description 1
- 230000035515 penetration Effects 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 238000000528 statistical test Methods 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
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- 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/58—Testing, adjusting or calibrating thereof
- A61B6/582—Calibration
- A61B6/583—Calibration using calibration phantoms
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/08—Projecting images onto non-planar surfaces, e.g. geodetic screens
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10081—Computed x-ray tomography [CT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20212—Image combination
- G06T2207/20216—Image averaging
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Medical Informatics (AREA)
- High Energy & Nuclear Physics (AREA)
- Molecular Biology (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Biophysics (AREA)
- Algebra (AREA)
- Operations Research (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Optics & Photonics (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- General Engineering & Computer Science (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Pulmonology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Probability & Statistics with Applications (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
Abstract
本公开实施例公开了一种基于钨靶X射线球管的双能CT系统的能谱数据校正方法，所述方法包括：利用能谱模拟软件模拟基于钨靶X射线球管发射的射线能谱，并作为初始能谱数据；基于初始能谱数据确定基于钨靶X射线球管发射的射线能谱的能谱拟合公式中的能谱参数，以及计算关键系数的迭代初值；基于初始能谱数据计算得到不同厚度衰减下的投影数据，作为校正初值；使用待校正的双能CT系统对模体组件进行扫描，实际采集得到不同厚度衰减下的投影数据，作为校正依据；基于校正初值、迭代初值和校正依据，使用迭代算法迭代修正关键系数；迭代结束后，基于关键系数更新能谱拟合公式；将能谱离散点带入所更新后的能谱拟合公式后，得到校正后的能谱数据。
Description
技术领域
本公开涉及X射线成像技术领域，具体涉及一种基于钨靶X射线球管的双能CT系统的能谱数据校正方法。
背景技术
CT成像的准确性依赖于对X射线与物质相互作用建模的准确性，其中X射线能谱数据的准确性是一个重要因素，如何获取X射线能谱信息一直是X射线CT成像中的一个重要研究课题。尤其在双能X射线CT成像技术中，双能CT成像技术可以实现物质识别和金属伪影去除，在临床医学上有广泛的应用和前景，而其重建过程中常用的基物质分解方法的主要原理是基于基物质(比如水和钙)在真实能谱下衰减情况的物理模型，因此，采用准确度更高的能谱数据进行衰减计算可以更符合真实衰减情况，进而实现更准确的基物质分解，最终提高双能CT成像技术中物质识别精度和图像质量。
在实际应用场景中，就算是技术参数相同的X射线球管(靶材料均为钨)，也会因为加工过程和使用环境无法做到一模一样，其电压波动、靶面不平滑程度和总滤过衰减效应无法完全一致，导致某单一的能谱数据不能适用于所有球管。
已有的X射线能谱获取方法可大致分为两类：直接方法和间接方法。直接方法即使用X射线光谱仪等设备直接测量X射线能谱；间接方法是利用X射线穿过不同厚度已知材料模体的测量数据与X射线能谱、模体厚度、模体衰减系数间的定量关系，建立关于X射线能谱的方程组，通过求解该方程组，获得X射线能谱信息。对于X射线CT系统，由于X光机产生的射线流强很大，利用光谱仪等设备直接测量X射线能谱在操作上有一定的难度，无法作为日常校正流程，而且能谱仪会增加成本。而间接测量方法中目前已有的方法主要是通过衰减测量实验来拟合能谱，在实际操作中，针对一种扫描条件(固定的电压、电流、滤波片等)，通常需要测量十几至几十个不同厚度模体的测量数据，工作繁复，所需时间长。
然而，间接测量方法中不同钨靶X射线球管发出的能谱曲线中，存在因电压波动、靶面不平滑、球管固有衰减情况不一致而带来的能谱曲线的偏移和波动，导致每台设备在双能CT成像中重建所用的能谱数据不一样的情况。
发明内容
本公开实施例提供一种基于钨靶X射线球管的双能CT系统的能谱数据校正方法。
第一方面，本公开实施例中提供了一种基于钨靶X射线球管的双能CT系统的能谱数据校正方法，包括：
利用能谱模拟软件模拟基于钨靶X射线球管发射的射线能谱，并作为初始能谱数据；
基于所述初始能谱数据确定所述基于钨靶X射线球管发射的射线能谱的能谱拟合公式中的4个能谱参数，以及计算所述能谱拟合公式中6个关键系数的迭代初值；
基于所述初始能谱数据计算得到不同厚度衰减下的投影数据，作为校正初值；
使用待校正的双能CT系统对模体组件进行扫描，实际采集得到所述不同厚度衰减下的投影数据，作为校正依据；
基于所述校正初值、6个关键系数的迭代初值和所述校正依据，使用迭代算法迭代修正所述6个关键系数；
迭代结束后，基于修正后的6个关键系数更新所述能谱拟合公式；
将能谱离散点带入所更新后的所述能谱拟合公式后，得到校正后的能谱数据；其中，所述能谱拟合公式如下表示：
其中，ΦkV(E)表示能谱数据，E为能谱数据中的能量离散点，c1、c2、c3、c4为4个能谱参数，根据初始能谱数据确定；a0、a1、a2、a3、bp、bq为6个关键系数；冲击函数δ(x·)用于描述钨靶的特征能量：57.98keV、59.32keV、67.2keV和69.1keV。
进一步地，基于所述初始能谱数据确定所述基于钨靶X射线球管发射的射线能谱的能谱拟合公式中的4个能谱参数，以及计算所述能谱拟合公式中6个关键系数的迭代初值，包括：
基于所述初始能谱数据确定能量在57.98keV、59.32keV、67.2keV和69.1keV下的数值，并基于所述能谱拟合公式以及所述能量在57.98keV、59.32keV、67.2keV和69.1keV下的数值求解得到4个能谱参数c1、c2、c3、c4；其中，系数bq的初始值为0；
将所述初始能谱数据中能量在57.98keV、59.32keV、67.2keV和69.1keV处的数值及前后邻近第一预定数量个点的数值去除，剩余数值进行插值得到没有特征峰的新能谱数据，且该新能谱数据满足曲线：
在能谱拟合公式中带入所述新能谱数据中的第二预定数量组能量离散点E和对应粒子数ΦkV(E)的数据，将能谱拟合公式中的5个关键系数a0、a1、a2、a3、bp求出，作为该5个关键系数a0、a1、a2、a3、bp系数的迭代初值。
进一步地，基于所述初始能谱数据计算得到不同厚度衰减下的投影数据，作为校正初值，包括：
利用如下公式基于所述初始能谱数据计算得到不同厚度衰减下的投影数据：
其中，
为投影数据，En为第n个能量离散点在初始能量数据中的能量值，ΦkV(En)为该能量值En下的粒子数，N为初始能谱数据的总离散点个数，μnl为第l个投影数据中所用的衰减材料在第n个能量离散点的衰减系数，Hl为第l个投影数据中所用的衰减材料的厚度，均为已知的数值。
进一步地，使用待校正的双能CT系统对模体组件进行扫描，实际采集得到所述不同厚度衰减下的投影数据，作为校正依据，包括：
将模体组件A和模体组件B分别放置在待校正的双能CT系统中进行扫描，获取两个第一图像；所述模体组件A和模体组件B上分别设置有多于5个的多个圆孔，所述多个圆孔内填充有不同厚度衰减系数的材质；
将所述模体组件A和模体组件B移除后，获取相同曝光条件下的第二图像；
基于所述第一图像中多个圆孔区域中的图像均值以及所述第二图像的典型区域均值获得所述不同厚度衰减下的投影数据，作为校正依据。
进一步地，基于所述校正初值、6个关键系数的迭代初值和所述校正依据，使用迭代算法迭代修正所述6个关键系数，包括：
将6个关键系数的迭代初值带入迭代算法构建的方程组中，所述方程组中的未知数为6个关键系数的修正值；
求解得到所述6个关键系数的修正值，并利用修正值更新得到所述6个关键系数的迭代初值；
利用更新后的所述6个关键系数的迭代初值以及能谱拟合公式更新所述校正初值；
基于更新后的所述6个关键系数的迭代初值、更新后的所述校正初值以及所述校正依据进行下一轮迭代，直到满足迭代终止条件为止。
进一步地，模体组件包括多个中空部分，每个中空部分填充有不同厚度的材料，所述中空部分横截面形状包括圆形、菱形、长条形中的一种或多种的组合；所述不同厚度的材料的厚度梯度基于待校正的双能CT系统的最大衰减厚度确定。
第二方面，本公开实施例中提供了一种基于钨靶X射线球管的双能CT系统的能谱数据校正装置，包括：
模拟模块，被配置为利用能谱模拟软件模拟基于钨靶X射线球管发射的射线能谱，并作为初始能谱数据；
确定模块，被配置为基于所述初始能谱数据确定所述基于钨靶X射线球管发射的射线能谱的能谱拟合公式中的4个能谱参数，以及计算所述能谱拟合公式中6个关键系数的迭代初值；
计算模块，被配置为基于所述初始能谱数据计算得到不同厚度衰减下的投影数据，作为校正初值；
扫描模块，被配置为使用待校正的双能CT系统对模体组件进行扫描，实际采集得到所述不同厚度衰减下的投影数据，作为校正依据；
修正模块，被配置为基于所述校正初值、6个关键系数的迭代初值和所述校正依据，使用迭代算法迭代修正所述6个关键系数；
更新模块，被配置为迭代结束后，基于修正后的6个关键系数更新所述能谱拟合公式；
获得模块，被配置为将能谱离散点带入所更新后的所述能谱拟合公式后，得到校正后的能谱数据；其中，所述能谱拟合公式如下表示：
其中，ΦkV(E)表示能谱数据，E为能谱数据中的能量离散点，c1、c2、c3、c4为4个能谱参数，根据初始能谱数据确定；a0、a1、a2、a3、bp、bq为6个关键系数；冲击函数δ(x·)用于描述钨靶的特征能量：57.98keV、59.32keV、67.2keV和69.1keV。
所述功能可以通过硬件实现，也可以通过硬件执行相应的软件实现。所述硬件或软件包括一个或多个与上述功能相对应的模块。
在一个可能的设计中，上述装置的结构中包括存储器和处理器，所述存储器用于存储一条或多条支持上述装置执行上述对应方法的计算机指令，所述处理器被配置为用于执行所述存储器中存储的计算机指令。上述装置还可以包括通信接口，用于上述装置与其他设备或通信网络通信。
第三方面，本公开实施例提供了一种电子设备，包括存储器、处理器以及存储在存储器上的计算机程序，其中，所述处理器执行所述计算机程序以实现上述任一方面所述的方法。
第四方面，本公开实施例提供了一种计算机可读存储介质，用于存储上述任一装置所用的计算机指令，该计算机指令被处理器执行时用于实现上述任一方面所述的方法。
第五方面，本公开实施例提供了一种计算机程序产品，其包含计算机指令，该计算机指令被处理器执行时用于实现上述任一方面所述的方法。
本公开实施例提供的技术方案可以包括以下有益效果：
本公开提出了一种双能CT系统专用能谱校正模体；基于该校正模体的扫描结果，本公开还提出了一种配套的校正流程。其中校正流程将能谱拟合公式表达为含有6个关键系数的公式(现有技术中为9个或7个)。在校正流程的算法中，因为非线性最小二乘算法迭代次数和结果对初值敏感，因此采用蒙特卡洛模拟的能谱数据来计算的投影值和能谱公式中6个关键系数值，将其作为后续修正迭代算法的初始值(现有技术中没有此做法，有论文使用开源能谱程序结果作为修正初始值，且计算过程不一样)。
应当理解的是，以上的一般描述和后文的细节描述仅是示例性和解释性的，并不能限制本公开。
附图说明
结合附图，通过以下非限制性实施方式的详细描述，本公开的其它特征、目的和优点将变得更加明显。在附图中：
图1示出根据本公开一实施方式的基于钨靶X射线球管的双能CT系统的能谱数据校正方法的流程图；
图2(a)和图2(b)示出根据本公开一实施方式的两种模体组件的侧面和横截面示意图；
图3示出根据本公开一实施方式的待校正双能CT系统的扫描结构示意图；
图4是适于用来实现根据本公开一实施方式的基于钨靶X射线球管的双能CT系统的能谱数据校正方法的电子设备的结构示意图。
具体实施方式
下文中，将参考附图详细描述本公开的示例性实施方式，以使本领域技术人员可容易地实现它们。此外，为了清楚起见，在附图中省略了与描述示例性实施方式无关的部分。
在本公开中，应理解，诸如“包括”或“具有”等的术语旨在指示本说明书中所公开的特征、数字、步骤、行为、部件、部分或其组合的存在，并且不排除一个或多个其他特征、数字、步骤、行为、部件、部分或其组合存在或被添加的可能性。
另外还需要说明的是，在不冲突的情况下，本公开中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本公开。
本公开提出了一种专用模体和配套校准算法，基于首先利用蒙特卡洛方法模拟出已知技术参数下球管的能谱，将其作为校准的初始能谱数据；再利用专用模体对特定双能CT设备进行投影测量，根据投影图中不同衰减厚度下的图像值，并将初始能谱数值作为迭代初值，结合迭代算法进行能谱函数的修正迭代，得到更新后的6个关键系数值。用于校准的能谱拟合公式是基于现技术领域已有的研究并结合X射线产生物理机制原理及特征进行了改进，将有效修正系数简化到6种，从而实现高效实用的能谱修正。
下面通过具体实施例详细介绍本公开实施例的细节。
图1示出根据本公开一实施方式的基于钨靶X射线球管的双能CT系统的能谱数据校正方法的流程图。如图1所示，该基于钨靶X射线球管的双能CT系统的能谱数据校正方法包括以下步骤：
在步骤S101中，利用能谱模拟软件模拟基于钨靶X射线球管发射的射线能谱，并作为初始能谱数据；
在步骤S102中，基于所述初始能谱数据确定所述基于钨靶X射线球管发射的射线能谱的能谱拟合公式中的4个能谱参数，以及计算所述能谱拟合公式中6个关键系数的迭代初值；
在步骤S103中，基于所述初始能谱数据计算得到不同厚度衰减下的投影数据，作为校正初值；
在步骤S104中，使用待校正的双能CT系统对模体组件进行扫描，实际采集得到所述不同厚度衰减下的投影数据，作为校正依据；
在步骤S105中，基于所述校正初值、6个关键系数的迭代初值和所述校正依据，使用迭代算法迭代修正所述6个关键系数；
在步骤S106中，迭代结束后，基于修正后的6个关键系数更新所述能谱拟合公式；
在步骤S107中，将能谱离散点带入所更新后的所述能谱拟合公式后，得到校正后的能谱数据；其中，所述能谱拟合公式如下表示：
其中，ΦkV(E)表示能谱数据，E为能谱数据中的能量离散点，c1、c2、c3、c4为4个能谱参数，根据初始能谱数据确定；a0、a1、a2、a3、bp、bq为6个关键系数；冲击函数δ(x·)用于描述钨靶的特征能量：57.98keV、59.32keV、67.2keV和69.1keV。
在本实施例的一个可选实现方式中，步骤S102，即基于所述初始能谱数据确定所述基于钨靶X射线球管发射的射线能谱的能谱拟合公式中的4个能谱参数，以及计算所述能谱拟合公式中6个关键系数的迭代初值的步骤，进一步包括以下步骤：
基于所述初始能谱数据确定能量在57.98keV、59.32keV、67.2keV和69.1keV下的数值，并基于所述能谱拟合公式以及所述能量在57.98keV、59.32keV、67.2keV和69.1keV下的数值求解得到4个能谱参数c1、c2、c3、c4；其中，系数bq的初始值为0；
将所述初始能谱数据中能量在57.98keV、59.32keV、67.2keV和69.1keV处的数值及前后邻近第一预定数量个点的数值去除，剩余数值进行插值得到没有特征峰的新能谱数据，且该新能谱数据满足曲线：
在能谱拟合公式中带入所述新能谱数据中的第二预定数量组能量离散点E和对应粒子数ΦkV(E)的数据，将能谱拟合公式中的5个关键系数a0、a1、a2、a3、bp求出，作为该5个关键系数a0、a1、a2、a3、bp系数的迭代初值。
在本实施例的一个可选实现方式中，步骤S103，即基于所述初始能谱数据计算得到不同厚度衰减下的投影数据，作为校正初值的步骤，进一步包括以下步骤：
利用如下公式基于所述初始能谱数据计算得到不同厚度衰减下的投影数据：
其中，
为投影数据，En为第n个能量离散点在初始能量数据中的能量值，ΦkV(En)为该能量值En下的粒子数，N为初始能谱数据的总离散点个数，μnl为第l个投影数据中所用的衰减材料在第n个能量离散点的衰减系数，Hl为第l个投影数据中所用的衰减材料的厚度，均为已知的数值。
在本实施例的一个可选实现方式中，步骤S104，即使用待校正的双能CT系统对模体组件进行扫描，实际采集得到所述不同厚度衰减下的投影数据，作为校正依据的步骤，进一步包括以下步骤：
将模体组件A和模体组件B分别放置在待校正的双能CT系统中进行扫描，获取两个第一图像；所述模体组件A和模体组件B上分别设置有多于5个的多个圆孔，所述多个圆孔内填充有不同厚度衰减系数的材质；
将所述模体组件A和模体组件B移除后，获取相同曝光条件下的第二图像；
基于所述第一图像中多个圆孔区域中的图像均值以及所述第二图像的典型区域均值获得所述不同厚度衰减下的投影数据，作为校正依据。
在本实施例的一个可选实现方式中，步骤S105，即基于所述校正初值、6个关键系数的迭代初值和所述校正依据，使用迭代算法迭代修正所述6个关键系数的步骤，进一步包括以下步骤：
将6个关键系数的迭代初值带入迭代算法构建的方程组中，所述方程组中的未知数为6个关键系数的修正值；
求解得到所述6个关键系数的修正值，并利用修正值更新得到所述6个关键系数的迭代初值；
利用更新后的所述6个关键系数的迭代初值以及能谱拟合公式更新所述校正初值；
基于更新后的所述6个关键系数的迭代初值、更新后的所述校正初值以及所述校正依据进行下一轮迭代，直到满足迭代终止条件为止。
下面通过一种具体实现方式详细介绍本公开的技术细节。
本公开首先对现有的基于钨靶X射线球管的能谱曲线拟合公式进行改进。X射线球管光谱分布与管电压、电流及靶材料的原子序数有关，相关研究提出的光谱模型由三部分构成：靶面产生的韧致辐射、特征辐射及球管固有滤过对射线的衰减。理论研究表明，靶面产生的韧致辐射强度会受各种因素影响，如电压波动、靶面不平滑等；球管固有滤过是指采用钨靶的球管产生的X射线在射出球管窗口时被球管玻壳、管芯与管套间的绝缘油、管套等衰减的总和，在医学诊断光谱段，物质对X线的衰减主要以光电效应和康普顿散射为主；第三部分的特征辐射则是与靶物质类型相关，对于钨材料，特征辐射发生在能量分别为57.98、59.32、67.20、69.10keV处，可用冲击函数表示。其能谱拟合公式如下：
其中，E为各能量离散点(假如双能CT的曝光电压值为120kV，则能谱中能量离散点可选取0～120keV中的120个等距分布点)，前面部分的三次多项式：(a0(ky)+a1(ky)E+a2(ky)E2+a3(ky)E3)用于表达轫致辐射能谱曲线；固有滤过的衰减用指数函数表达：
后面的4个冲击函数δ(x)用于描述钨靶的特征能量：57.98keV、59.32keV、67.2keV和69.1keV。上式为目前该技术领域常用的能谱曲线拟合公式，共有9个参数a0、a1、a2、a3、bp、c1、c2、c3、c4。在特征辐射表达部分，对于都是基于钨靶X射线球管的双能CT系统来说，特征辐射能量(57.98keV、59.32keV、67.2keV和69.1keV)是一致的，但其强度与球管固有滤过有关，因此对上式进行改进，引入特征能量强度根据固有滤过衰减效应改变的参数bq，拟合公式改进为：
其中，参数c1、c2、c3、c4将根据初始能谱曲线确定，校准流程中修正的参数为a0、a1、a2、a3、bp、bq共6个系数；冲击函数δ(·)用于描述钨靶的特征能量：57.98keV、59.32keV、67.2keV和69.1keV。
对于不同的球管，因为加工过程和使用环境无法做到一模一样，其电压波动、靶面不平滑程度和总滤过衰减效应是无法达到完全一致的，这些因素也最终导致某单一的能谱数据不能适用于所有采用钨靶的X射线球管。因此，基于上式，引入不同球管带来的不稳定因素，对涉及变动的6个参数a0、a1、a2、a3、bp、bq用实际采集数据进行修正。
为了修正上述能谱公式，本公开提供了一组专用模体，包括两块圆柱形的模体组件(包括A组件和B组件)，每块中间都有10个小的中空圆柱，A组件中的10个中空圆柱由不同厚度的铜材料进行填充，B组件中的10个中空圆柱是由不同厚度的铝材料进行填充(如图1、图2所示)。结合临床用双能CT设备的电压值分布在80kV～140kV之间，模体组件中铜材料和铝材料选用的厚度梯度是基于各曝光电压下最大衰减厚度(大于这个厚度射线穿透较少影响测量)，确保每个电压下的模体组件拍摄图像中都至少有5个小圆孔里的衰减图像值是可用的。需要说明的是，双能CT是一种指用两个高、低电压切换曝光的CT，上文中的各曝光电压指的双能CT的高能曝光电压和低能曝光电压。因此，在确定了曝光电压，比如高能曝光电压下(比如140kV)铜物质的最大衰减厚度是5mm，那么圆孔厚度在小于5mm的范围内的圆孔里的衰减图像值是可用的，所以设计圆孔厚度时确保小于5mm厚度范围内至少有5个孔即可。
其中，最大衰减厚度的计算公式如下：
其中，Hmax为最大衰减厚度；Emax为X射线光子的最大能量，由曝光电压值确定；μ(Emax/2)为圆孔中的物质在Emax/2这个能量下的线性衰减系数，均为业内查表可知的已知值。
将蒙特卡洛模拟所得能谱数据计算得到的与实验中相同厚度衰减下的投影数据作为修正的初值，再利用真实采集的专用模体中不同厚度衰减下的投影数据作为校正依据，最后采用非线性最小二乘迭代算法进行计算来修正能谱公式中的待求参数。需要说明的是，迭代算法可以采用除非线性最小二乘之外的其他迭代算法，具体可以根据实际应用进行选择，在此不做限制。
需要说明的是，本公开实施例中模体组件中圆孔的排列方式、间距等不限于图中所示，可以以任意方式进行排列，只要任意两孔之间不重叠即可。另外，模体组件的圆孔中的材料不限于铜和铝，可以替换成任意线衰减系数确定的物质。在一些实施例中，模体组件的圆孔中的材料为单种物质，该材料可以是金属材料和非金属材料，非金属材料例如可以是均匀的有机玻璃材质。
下面通过具体的实例介绍本公开的一种实现方式。
①蒙特卡洛模拟得到初始能谱数据
蒙特卡洛模拟方法，也称为计算机随机模拟方法、统计模拟法、统计试验法，是基于“随机数”的计算方法，或者是说把概率现象作为研究对象的数值模拟方法，目前广泛用于核物理、医学物理等领域。本公开提出的方法中先利用蒙特卡洛软件模拟双能CT系统所用技术参数下的基于钨靶的X射线球管的作用过程。首先在蒙特卡洛专用程序中建立与真实光机相同的物理模型(包括尺寸、衰减厚度、靶物质类型、大小、反应面积、电子能量、电子数量、出束口大小)，然后设定模拟的电子数量，该程序就会自动计算这么多数量的电子一个一个与靶物质反应产生的X射线光子的数量和能量，这些不同光子的能量组成的曲线就是能谱曲线(一般横坐标为能量，纵坐标为个数或归一化后的个数，也叫发生概率)；
模拟中建立物理模型的条件(电子能量)需与双能CT系统实际所用的高能和低能曝光电压参数(比如140kV和80kV)一致，且其中固有滤过在技术参数所给的数值减少1mmAl的厚度，用于预留后续修正空间。
上述操作下得到能谱数据，找到能量在57.98keV、59.32keV、67.2keV和69.1keV下的数值求得拟合能谱公式中的c1、c2、c3、c4。而系数bq的初始值为0，即
②获取能谱公式中关键系数初始值
根据①中的模拟数据曲线，去掉特征能量下57.98keV、59.32keV、67.2keV和69.1keV处及附近前后5个点的数值，将剩下的数值重新插值(公式如下(3)所示)得到没有特征峰的新能谱数据，且其满足曲线：
因此在数据中代入几组(大于5组即可)能量E和对应粒子数ΦkV(E)的数据，将公式中的系数a0、a1、a2、a3、bp求出，作为系数初始值。
上面提到的插值公式如下：
(x,y)为插值点在能谱曲线中的坐标，(x1,y1)、(x2,y2)为已有的前后近邻两点在能谱曲线中的点。
③获取模体投影数据
用待校正的双能CT系统对专用模体进行扫描，其中曝光电压与蒙特卡洛模拟中一致，即为该系统正常使用中所用高能和低能曝光电压(比如140kV和80kV)。具体操作为分别将模体组件A和B(如图2(a)-图2(b)所示)放在X射线发生器和探测器之间，确保探测器影像能涵盖模体的所有部分(如图3所示)，并获取第一图像；然后再移除模体组件，获取同样曝光条件下的第二图像；在模体组件A和B的第一图像中各选择5个小圆孔区域作为统计区域，模体投影数据计算如下：
Pl＝-ln(Il/I0)，l＝1，2，3....10 (4)
Il为第一图像中第l个统计区域中的图像均值，I0为去掉模体扫描后的第二图像中图像典型区域均值。
从上述公式可以看出，从模体组件A和B分别获得5个圆孔区域的图像均值，总共是10个统计区域的图像均值，而实际应用中，可以分别获得多于5个圆孔区域的图像均值，也即总共可以获得多于10个统计区域的图像均值。
④迭代修正计算
该流程以高能或低能中任意一个能谱的修正为例，另一个的修正过程与其完全一致，只是基于其对应的扫描数据即可。
为了方便，将能谱公式中的待求参数a0、a1、a2、a3、bp、bq记为xm,m＝1,2,..6；再将实际采集得到的投影数据记为Pl,l＝1,2,3....10，将利用能谱数据计算得到的投影数据记为
计算方式如下：
其中，En为第n个能量点的能量值，ΦkV(Fn)为该能量下的粒子数，N为能谱数据的总离散点个数，这些数据通过蒙特卡洛计算的初始能谱已经获得。μnl为第l个投影数据中所用的衰减材料(铜或铝)在第n个能量点下的衰减系数，Hl为第l个投影数据中所用的衰减材料的厚度，均为已知的数值(当作常数项)，因此
还是对应6个未知数x1～x6的一个方程组。
迭代修正的思路就是基于真实测量得到的投影数据Pl和初始能谱数据计算下的投影数据
之差来修正能谱公式中的系数值，迭代算法步骤如下：
a)计算迭代初始值
代入蒙特卡洛模拟能谱结果和实验中的衰减厚度和各能量下的衰减系数信息到公式(5)中，计算
b)代入特定的迭代格式，以基于非线性最小二乘原理的Gauss-Newton迭代算法为例，其构造方程组：
其中：
代入6个待求参数值x1～x6的初始值(也即6个关键系数a0、a1、a2、a3、bp、bq，步骤①②中已求得)，方程组中除了Δ1～Δ6其他均为常数，因此直接求解上述方程组中的Δ1～Δ6；
c)设置迭代条件
上述方程组求解得到的Δl作为参数xl的修正值，即
其中K为迭代次数。每次求解得到一组Δl，用其修正原来的参数值xi，再将修正后的xi代入方程组(
和
也在修正后的xl的基础上更新一遍)再次求解Δi，如此反复循环。每次求解之后加入判别条件：
满足上述条件则迭代终止，否则重复步骤b)～c)。
⑤获取修正后的能谱
迭代计算结束后，得到最终的能谱拟合公式中的关键系数x1～x6，对应公式(2)中的参数a0、a1、a2、a3、bp、bq，代入能谱离散点：比如双能CT系统高能下的扫描条件为140kV，则选取0～140keV能量的200个等距能量值代入能谱公式进行计算；由此获得修正后的能谱数据。
下述为本公开装置实施例，可以用于执行本公开方法实施例。
根据本公开一实施方式的基于钨靶X射线球管的双能CT系统的能谱数据校正装置，该装置可以通过软件、硬件或者两者的结合实现成为电子设备的部分或者全部。该基于钨靶X射线球管的双能CT系统的能谱数据校正装置包括：
模拟模块，被配置为利用能谱模拟软件模拟基于钨靶X射线球管发射的射线能谱，并作为初始能谱数据；
确定模块，被配置为基于所述初始能谱数据确定所述基于钨靶X射线球管发射的射线能谱的能谱拟合公式中的4个能谱参数，以及计算所述能谱拟合公式中6个关键系数的迭代初值；
计算模块，被配置为基于所述初始能谱数据计算得到不同厚度衰减下的投影数据，作为校正初值；
扫描模块，被配置为使用待校正的双能CT系统对模体组件进行扫描，实际采集得到所述不同厚度衰减下的投影数据，作为校正依据；
修正模块，被配置为基于所述校正初值、6个关键系数的迭代初值和所述校正依据，使用迭代算法迭代修正所述6个关键系数；
更新模块，被配置为迭代结束后，基于修正后的6个关键系数更新所述能谱拟合公式；
获得模块，被配置为将能谱离散点带入所更新后的所述能谱拟合公式后，得到校正后的能谱数据；其中，所述能谱拟合公式如下表示：
其中，ΦkV(E)表示能谱数据，E为能谱数据中的能量离散点，c1、c2、c3、c4为4个能谱参数，根据初始能谱数据确定；a0、a1、a2、a3、bp、bq为6个关键系数；冲击函数δ(x·)用于描述钨靶的特征能量：57.98keV、59.32keV、67.2keV和69.1keV。
本实施例中的基于钨靶X射线球管的双能CT系统的能谱数据校正装置与上文中的基于钨靶X射线球管的双能CT系统的能谱数据校正方法对应一致，具体细节可以参见上文中对基于钨靶X射线球管的双能CT系统的能谱数据校正方法的描述，在此不再赘述。图4是适于用来实现根据本公开一实施方式的基于钨靶X射线球管的双能CT系统的能谱数据校正方法的电子设备的结构示意图。
如图4所示，电子设备400包括处理单元401，其可实现为CPU、GPU、FPGA、NPU等处理单元。处理单元401可以根据存储在只读存储器(ROM)402中的程序或者从存储部分408加载到随机访问存储器(RAM)403中的程序而执行本公开上述任一方法的实施方式中的各种处理。在RAM403中，还存储有电子设备400操作所需的各种程序和数据。处理单元401、ROM402以及RAM403通过总线404彼此相连。输入/输出(I/O)接口405也连接至总线404。
以下部件连接至I/O接口405：包括键盘、鼠标等的输入部分406；包括诸如阴极射线管(CRT)、液晶显示器(LCD)等以及扬声器等的输出部分407；包括硬盘等的存储部分408；以及包括诸如LAN卡、调制解调器等的网络接口卡的通信部分409。通信部分409经由诸如因特网的网络执行通信处理。驱动器410也根据需要连接至I/O接口405。可拆卸介质411，诸如磁盘、光盘、磁光盘、半导体存储器等等，根据需要安装在驱动器410上，以便于从其上读出的计算机程序根据需要被安装入存储部分408。
特别地，根据本公开的实施方式，上文参考本公开实施方式中的任一方法可以被实现为计算机软件程序。例如，本公开的实施方式包括一种计算机程序产品，其包括有形地包含在机器可读介质上的计算机程序，所述计算机程序包含用于执行本公开实施方式中任一方法的程序代码。在这样的实施方式中，该计算机程序可以通过通信部分409从网络上被下载和安装，和/或从可拆卸介质411被安装。
附图中的流程图和框图，图示了按照本公开各种实施方式的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，路程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个接连地表示的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或操作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
描述于本公开实施方式中所涉及到的单元或模块可以通过软件的方式实现，也可以通过硬件的方式来实现。所描述的单元或模块也可以设置在处理器中，这些单元或模块的名称在某种情况下并不构成对该单元或模块本身的限定。
作为另一方面，本公开还提供了一种计算机可读存储介质，该计算机可读存储介质可以是上述实施方式中所述装置中所包含的计算机可读存储介质；也可以是单独存在，未装配入设备中的计算机可读存储介质。计算机可读存储介质存储有一个或者一个以上程序，所述程序被一个或者一个以上的处理器用来执行描述于本公开的方法。
以上描述仅为本公开的较佳实施例以及对所运用技术原理的说明。本领域技术人员应当理解，本公开中所涉及的发明范围，并不限于上述技术特征的特定组合而成的技术方案，同时也应涵盖在不脱离所述发明构思的情况下，由上述技术特征或其等同特征进行任意组合而形成的其它技术方案。例如上述特征与本公开中公开的(但不限于)具有类似功能的技术特征进行互相替换而形成的技术方案。
Claims (7)
1.一种基于钨靶X射线球管的双能CT系统的能谱数据校正方法，其特征在于，包括：
利用能谱模拟软件模拟基于钨靶X射线球管发射的射线能谱，并作为初始能谱数据；
基于所述初始能谱数据确定所述基于钨靶X射线球管发射的射线能谱的能谱拟合公式中的4个能谱参数，以及计算所述能谱拟合公式中6个关键系数的迭代初值；
基于所述初始能谱数据计算得到不同厚度衰减下的投影数据，作为校正初值；
使用待校正的双能CT系统对模体组件进行扫描，实际采集得到所述不同厚度衰减下的投影数据，作为校正依据；
基于所述校正初值、6个关键系数的迭代初值和所述校正依据，使用迭代算法迭代修正所述6个关键系数；
迭代结束后，基于修正后的6个关键系数更新所述能谱拟合公式；
将能谱离散点带入更新后的所述能谱拟合公式后，得到校正后的能谱数据；其中，所述能谱拟合公式如下表示：
其中，ΦkV(E)表示能谱数据中的粒子数，E为能谱数据中的能量离散点，c1、c2、c3、c4为4个能谱参数，根据初始能谱数据确定；a0、a1、a2、a3、bp、bq为6个关键系数；冲击函数δ(x·)用于描述钨靶的特征能量：57.98keV、59.32keV、67.2keV和69.1keV；
基于所述初始能谱数据确定所述基于钨靶X射线球管发射的射线能谱的能谱拟合公式中的4个能谱参数，以及计算所述能谱拟合公式中6个关键系数的迭代初值，包括：
基于所述初始能谱数据确定能量在57.98keV、59.32keV、67.2keV和69.1keV下的数值，
并基于所述能谱拟合公式以及所述能量在57.98keV、59.32keV、67.2keV和69.1keV下的数值求解得到4个能谱参数c1、c2、c3、c4；其中，关键系数bq的初始值为0；
将所述初始能谱数据中能量在57.98keV、59.32keV、67.2keV和69.1keV处的数值及前后邻近第一预定数量个点的数值去除，剩余数值进行插值得到没有特征峰的新能谱数据，且该新能谱数据满足曲线：
在能谱拟合公式中带入所述新能谱数据中的第二预定数量组能量离散点和对应粒子数的数据，将能谱拟合公式中的5个关键系数a0、a1、a2、a3、bp求出，作为该5个关键系数a0、a1、a2、a3、bp系数的迭代初值；
基于所述初始能谱数据计算得到不同厚度衰减下的投影数据，作为校正初值，包括：
利用如下公式基于所述初始能谱数据计算得到不同厚度衰减下的投影数据：
其中，
为投影数据；En为第n个能量离散点在初始能量数据中的能量值，ΦkV(En)为该能量值En下的粒子数，N为初始能谱数据的总离散点个数，μnl为第l个投影数据中所用的衰减材料在第n个能量离散点的衰减系数，Hl为第l个投影数据中所用的衰减材料的厚度，均为已知的数值。
2.根据权利要求1所述的方法，其特征在于，使用待校正的双能CT系统对模体组件进行扫描，实际采集得到所述不同厚度衰减下的投影数据，作为校正依据，包括：
将模体组件A和模体组件B分别放置在待校正的双能CT系统中进行扫描，获取两个第一图像；所述模体组件A和模体组件B上分别设置有多于5个的多个圆孔，所述多个圆孔内填充有不同厚度衰减系数的材质；
将所述模体组件A和模体组件B移除后，获取相同曝光条件下的第二图像；
基于所述第一图像中多个圆孔区域中的图像均值以及所述第二图像的典型区域均值获得所述不同厚度衰减下的投影数据，作为校正依据。
3.根据权利要求1所述的方法，其特征在于，基于所述校正初值、6个关键系数的迭代初值和所述校正依据，使用迭代算法迭代修正所述6个关键系数，包括：
将6个关键系数的迭代初值带入迭代算法构建的方程组中，所述方程组中的未知数为6个关键系数的修正值；
求解得到所述6个关键系数的修正值，并利用修正值更新得到所述6个关键系数的迭代初值；
利用更新后的所述6个关键系数的迭代初值以及能谱拟合公式更新所述校正初值；
基于更新后的所述6个关键系数的迭代初值、更新后的所述校正初值以及所述校正依据进行下一轮迭代，直到满足迭代终止条件为止。
4.根据权利要求1所述的方法，其特征在于，模体组件包括多个中空部分，每个中空部分填充有不同厚度的材料，所述中空部分横截面形状包括圆形、菱形、长条形中的一种或多种的组合；所述不同厚度的材料的厚度梯度基于待校正的双能CT系统的最大衰减厚度确定。
5.一种基于钨靶X射线球管的双能CT系统的能谱数据校正装置，其特征在于，包括：
模拟模块，被配置为利用能谱模拟软件模拟基于钨靶X射线球管发射的射线能谱，并作为初始能谱数据；
确定模块，被配置为基于所述初始能谱数据确定所述基于钨靶X射线球管发射的射线能谱的能谱拟合公式中的4个能谱参数，以及计算所述能谱拟合公式中6个关键系数的迭代初值；
计算模块，被配置为基于所述初始能谱数据计算得到不同厚度衰减下的投影数据，作为校正初值；
扫描模块，被配置为使用待校正的双能CT系统对模体组件进行扫描，实际采集得到所述不同厚度衰减下的投影数据，作为校正依据；
修正模块，被配置为基于所述校正初值、6个关键系数的迭代初值和所述校正依据，使用迭代算法迭代修正所述6个关键系数；
更新模块，被配置为迭代结束后，基于修正后的6个关键系数更新所述能谱拟合公式；
获得模块，被配置为将能谱离散点带入更新后的所述能谱拟合公式后，得到校正后的能谱数据；其中，所述能谱拟合公式如下表示：
其中，ΦkV(E)表示能谱数据中的粒子数，E为能谱数据中的能量离散点，c1、c2、c3、c4为4个能谱参数，根据初始能谱数据确定；a0、a1、a2、a3、bp、bq为6个关键系数；冲击函数δ(x·)用于描述钨靶的特征能量：57.98keV、59.32keV、67.2keV和69.1keV；
其中，所述确定模块，具体包括：
基于所述初始能谱数据确定能量在57.98keV、59.32keV、67.2keV和69.1keV下的数值，并基于所述能谱拟合公式以及所述能量在57.98keV、59.32keV、67.2keV和69.1keV下的数值求解得到4个能谱参数c1、c2、c3、c4；其中，关键系数bq的初始值为0；
将所述初始能谱数据中能量在57.98keV、59.32keV、67.2keV和69.1keV处的数值及前后邻近第一预定数量个点的数值去除，剩余数值进行插值得到没有特征峰的新能谱数据，且该新能谱数据满足曲线：
在能谱拟合公式中带入所述新能谱数据中的第二预定数量组能量离散点和对应粒子数的数据，将能谱拟合公式中的5个关键系数a0、a1、a2、a3、bp求出，作为该5个关键系数a0、a1、a2、a3、bp系数的迭代初值；
所述计算模块，具体包括：
利用如下公式基于所述初始能谱数据计算得到不同厚度衰减下的投影数据：
其中，
为投影数据；En为第n个能量离散点在初始能量数据中的能量值，ΦkV(En)为该能量值En下的粒子数，N为初始能谱数据的总离散点个数，μnl为第l个投影数据中所用的衰减材料在第n个能量离散点的衰减系数，Hl为第l个投影数据中所用的衰减材料的厚度，均为已知的数值。
6.一种电子设备，其特征在于，包括存储器、处理器以及存储在存储器上的计算机程序，其特征在于，所述处理器执行所述计算机程序以实现权利要求1-4任一项所述的方法。
7.一种计算机可读存储介质，其上存储有计算机指令，其特征在于，该计算机指令被处理器执行时实现权利要求1-4任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210927488.7A CN115270075B (zh) | 2022-08-03 | 2022-08-03 | 基于钨靶x射线球管的双能ct系统的能谱数据校正方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210927488.7A CN115270075B (zh) | 2022-08-03 | 2022-08-03 | 基于钨靶x射线球管的双能ct系统的能谱数据校正方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115270075A CN115270075A (zh) | 2022-11-01 |
| CN115270075B true CN115270075B (zh) | 2023-04-25 |
Family
ID=83749947
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210927488.7A Active CN115270075B (zh) | 2022-08-03 | 2022-08-03 | 基于钨靶x射线球管的双能ct系统的能谱数据校正方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115270075B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117270024B (zh) * | 2023-11-20 | 2024-02-20 | 合肥综合性国家科学中心人工智能研究院(安徽省人工智能实验室) | 能谱响应函数的校正方法、装置、计算机设备及存储介质 |
| CN117372287B (zh) * | 2023-12-06 | 2024-03-15 | 赛诺威盛科技(北京)股份有限公司 | 探测器能谱响应的补偿校正、装置、设备和介质 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110811660A (zh) * | 2019-10-25 | 2020-02-21 | 赛诺威盛科技（北京）有限公司 | 一种校正ct射线束硬化伪影的方法 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10869646B2 (en) * | 2018-04-12 | 2020-12-22 | Canon Medical Systems Corporation | Method and apparatus for computed tomography (CT) and material decomposition with pile-up correction calibrated using real pulse pileup effect and detector response |
| CN111643104B (zh) * | 2020-02-28 | 2022-03-18 | 清华大学 | Ct散射校正方法及系统 |
| CN111968060B (zh) * | 2020-08-28 | 2022-07-08 | 首都师范大学 | 一种基于倾斜投影修正技术的多能谱ct快速迭代重建方法 |
-
2022
- 2022-08-03 CN CN202210927488.7A patent/CN115270075B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110811660A (zh) * | 2019-10-25 | 2020-02-21 | 赛诺威盛科技（北京）有限公司 | 一种校正ct射线束硬化伪影的方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115270075A (zh) | 2022-11-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN115270075B (zh) | 基于钨靶x射线球管的双能ct系统的能谱数据校正方法 | |
| Wang et al. | Acuros CTS: A fast, linear Boltzmann transport equation solver for computed tomography scatter–Part II: System modeling, scatter correction, and optimization | |
| JP5264268B2 (ja) | 単位面積質量画像の作成方法 | |
| CN107802280B (zh) | 校正曲线生成方法、投影图像的校正方法、系统及存储介质 | |
| Ducote et al. | Scatter correction in digital mammography based on image deconvolution | |
| CN110175957B (zh) | 多能量ct基材料物质分解方法 | |
| JP2008502395A (ja) | 投影放射線撮影、特に乳房撮影における散乱放射線補正装置および方法 | |
| Faulkner et al. | Analysis of x-ray computed tomography images using the noise power spectrum and autocorrelation function | |
| Ding et al. | Image‐domain multimaterial decomposition for dual‐energy CT based on prior information of material images | |
| US10605933B2 (en) | X-ray spectral calibration technique for cone-beam CT | |
| Starman et al. | Investigation into the optimal linear time‐invariant lag correction for radar artifact removal | |
| Zhao et al. | Robust beam hardening artifacts reduction for computed tomography using spectrum modeling | |
| US20180328865A1 (en) | Apparatus for determining an effective energy spectrum of an x-ray tube | |
| Hayes et al. | Low‐dose cone‐beam CT via raw counts domain low‐signal correction schemes: Performance assessment and task‐based parameter optimization (Part I: Assessment of spatial resolution and noise performance) | |
| CN109916933B (zh) | 基于卷积神经网络的x射线计算机断层成像能谱估计方法 | |
| Marshall et al. | Performance evaluation of digital breast tomosynthesis systems: physical methods and experimental data | |
| Kontson et al. | Bowtie filters for dedicated breast CT: theory and computational implementation | |
| JP2017221339A (ja) | Ｘ線ｃｔ画像再構成方法およびコンピュータプログラム | |
| Talbi et al. | Local radiation dosimetry method using optically stimulated pulsed luminescence and Monte Carlo simulation | |
| Bradley | Estimation of bias and variance of measurements made from tomography scans | |
| Yu et al. | Heel effect adaptive flat field correction of digital x‐ray detectors | |
| van der Heyden et al. | Virtual monoenergetic micro-CT imaging in mice with artificial intelligence | |
| Leon et al. | Characterization of scatter in digital mammography from use of Monte Carlo simulations and comparison to physical measurements | |
| Shamshad | Characterisation of the x-ray source and camera in the MARS spectral system | |
| JP2017113526A (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 |