CN117665900A - 一种x射线图像校正方法及系统 - Google Patents
一种x射线图像校正方法及系统 Download PDFInfo
- Publication number
- CN117665900A CN117665900A CN202311682374.1A CN202311682374A CN117665900A CN 117665900 A CN117665900 A CN 117665900A CN 202311682374 A CN202311682374 A CN 202311682374A CN 117665900 A CN117665900 A CN 117665900A
- Authority
- CN
- China
- Prior art keywords
- energy
- ray
- projection value
- geometric body
- photon
- 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.)
- Granted
Links
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
-
- 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
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
- G01T1/00—Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
- G01T1/36—Measuring spectral distribution of X-rays or of nuclear radiation spectrometry
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
- G01T7/00—Details of radiation-measuring instruments
- G01T7/005—Details of radiation-measuring instruments calibration techniques
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N7/00—Computing arrangements based on specific mathematical models
- G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/08—Probabilistic or stochastic CAD
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; 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/10116—X-ray image
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Molecular Biology (AREA)
- High Energy & Nuclear Physics (AREA)
- Medical Informatics (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Spectroscopy & Molecular Physics (AREA)
- General Engineering & Computer Science (AREA)
- Optics & Photonics (AREA)
- Surgery (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Veterinary Medicine (AREA)
- Public Health (AREA)
- General Health & Medical Sciences (AREA)
- Animal Behavior & Ethology (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Computer Hardware Design (AREA)
- Biophysics (AREA)
- Heart & Thoracic Surgery (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Biomedical Technology (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Probability & Statistics with Applications (AREA)
- Algebra (AREA)
- Artificial Intelligence (AREA)
- Software Systems (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Computing Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Measurement Of Radiation (AREA)
Abstract
本发明公开了一种X射线图像校正方法及系统,所述X射线图像校正方法包括:模拟X射线球管管电压能量为E的X射线能谱分布;按照能谱划分步长将X射线能谱分布划分成若干个单能的子能谱;模拟各个子能谱X射线穿过不同贯穿长度L的被检测物体的光子通量;将子能谱的光子通量转换成投影值P;将贯穿长度L与各子能谱加权后的投影值P进行多项式拟合;利用得到贯穿长度L和投影值P的拟合曲线校正实际的多能投影值Pm,得到理想的单能投影值Ps,实现X射线图像的硬化校正。本发明的一种X射线图像校正方法,基于蒙卡模拟克服了硬化校正对实验条件所限制的问题,具有灵活性高、模拟精度高的优势,能有效提高重建图像的均匀性。
Description
技术领域
本发明涉及医学影像硬化校正技术领域,具体的涉及一种X射线图像校正方法及系统。
背景技术
对于理想的单能射线源,在穿过均匀物体时,射线衰减遵守Beer定理,有
I=I0e-μL
变换得
其中,I表示穿透物体后的射线强度,I0表示入射射线强度,μ表示被测物体对于射线的线性衰减系数,L表示射线透射物体的贯穿长度,Ps称为单能投影值。在单能情况下单能投影值与射线经过被检测物体的厚度成线性关系。
多能情形时,投影值表示为
其中,I′表示穿透物体后的射线强度,表示入射射线平均强度,μEi表示被测物体在射线能量为Ei的线性衰减系数。当X射线具有连续能谱分布时,由于物体的衰减,X射线能谱变硬,其穿透能力增强,投影值与射线贯穿被检测物体的长度成非线性关系,多能X射线的投影值小于单能时的投影值,且随着贯穿长度的增加,与单能情形下的投影值相差越大。图1为单能和多能时,投影值与贯穿长度的函数曲线。
因此,硬化校正目的在于将多能群X射线的投影值Pm校正为单能情形时的投影值Ps,通过建立投影值Pm与L之间的函数关系,通过中间量贯穿长度L,可计算单能的理想情形的投影值Ps,作为CT或CBCT图像重建输入。
传统的射束硬化校正方法,需要针对每一种材料试验测量该材料对射线的投影值与贯穿长度曲线。因此投影值与贯穿长度曲线对实验条件有很大的依赖性,每当改变X光机电压或者被测工件的材料等条件时,需要重新测量曲线才能完成硬化校正过程,这种方法费事费时。
为了突破硬化校正对实验条件依赖的限制,特提出本发明专利。
发明内容
为了解决上述技术问题,本发明提供一种X射线图像校正方法、系统及存储介质,X射线图像校正方法根据能谱划分后各子能谱加权后的投影值与贯穿被检测物体长度进行多项式拟合硬化校正,具体地,采用了如下技术方案:
一种X射线图像校正方法,包括:
模拟X射线球管管电压能量为E的X射线能谱分布;
按照能谱划分步长将X射线能谱分布划分成若干个单能的子能谱;
模拟各个子能谱X射线穿过不同贯穿长度L的被检测物体的光子通量;
将子能谱的光子通量转换成投影值P;
将贯穿长度L与各子能谱加权后的投影值P进行多项式拟合;
利用得到贯穿长度L和投影值P的拟合曲线校正实际的多能投影值Pm,得到理想的单能投影值Ps,实现X射线图像的硬化校正。
作为本发明的可选实施方式,本发明的一种X射线图像校正方法中,所述模拟各个子能谱X射线穿过不同贯穿长度L的被检测物体的光子通量利用蒙卡模拟,包括:
光子通量分布模拟几何建模,创建表示光子源的第一几何体,创建表示被检测物体的第二几何体,创建表示探测器平板的第三几何体和第四几何体,所述的第二几何体、第三几何体和第四几何体依次排布在第一几何体的一侧,所述第四几何体布置在距离第一几何体Smax位置处,其中,所述Smax表征距离光子源无穷远位置,所述第二几何体的厚度设置表征被检测物体的贯穿长度L;
光子通量分布模拟物理建模,为第二几何体、第三几何体和第四几何体模型选择相应的材料,源建模设置源粒子类型为光子源,光子源为固定点源并设置点的坐标值,设置光子源的能量为单能,大小为X射线能谱分布划分成若干个单能的子能谱的大小,设置光子源的方向参数为各向异性,光子源单方向出射;
光子通量分布模拟计数建模,计数粒子类型为光子,计数栅元号为用于计数的几何模型栅元号,统计光子源单能射束穿过物体光子数。
作为本发明的可选实施方式,本发明的一种X射线图像校正方法,将利用蒙卡模拟的各个子能谱X射线穿过不同贯穿长度L的被检测物体的光子通量分别保存为单能通量,构建含有各单能通量的通量数据库,在进行X射线图像的硬化校正时,根据划分成的若干个单能的子能谱调用通量数据库中相应的单能通量,快速计算拟合曲线。
作为本发明的可选实施方式,本发明的一种X射线图像校正方法中,所述将子能谱的光子通量转换成投影值P包括:
将蒙卡模拟各单能X射线穿过被检测物体后得到光子通量,将单能射束模拟的光子通量转成探测器平板响应,利用探测器平板响应正比于射线强度,累加归一化后计算得到投影值,探测器平板响应计算公式为
其中,i表示能群,t表示厚度,μi表示第i个能群对应被检测物体的质量衰减系数,μen,i表示第i个能群对应CsI的质量能量吸收系数,Φi,0表示蒙卡模拟的第i个能群透过厚度为t的光子通量,表示第i个能群的平均能量,Ii,t表示第i个能群透过厚度为t的水层在CsI探测器吸收的能量,Ii,t正比于探测器平板响应,将i个能群的Ii,t进行累加做归一化处理得到It,表示加权后的多能X射束透过厚度为t的被检测物体在CsI测探测器吸收的能量,计算公式如下:
利用平板亮度计算投影值,公式如下:
作为本发明的可选实施方式,本发明的一种X射线图像校正方法中,所述将贯穿长度L与各子能谱加权后的投影值P进行多项式拟合包括:
采用投影值P作为自变量,材料的贯穿长度L作为因变量,进行N次多项式最小二乘法拟合,得到L与P的函数关系如下:
其中,N不小于4。
作为本发明的可选实施方式,本发明的一种X射线图像校正方法中,所述利用得到贯穿长度L和投影值P的拟合曲线校正实际的多能投影值Pm,得到理想的单能投影值Ps,采用如下公式计算:
Ps=kf-1(Pm)
作为CT或CBCT等图像的多项式硬化校正输入。
作为本发明的可选实施方式,本发明的一种X射线图像校正方法中,所述模拟X射线球管管电压能量为E的X射线能谱分布利用蒙卡模拟,包括:
几何建模,创建表示电子源的第五几何体,创建表示钨靶的第六几何体,创建表示铝滤过板的第七几何体,创建表示探测器平板的第八几何体,所述的第六几何体具有靶角θ的倾斜靶面,所述的第五几何体位于第六几何体的倾斜靶面的竖直上侧,所述的第七几何体、第八几何体依次布置在第六几何体的倾斜靶面的水平右侧;
物理建模,针对本模拟实验场景使用的第六几何体、第七几何体、第八几何体选择相应的材料进行材料建模,源建模设置源粒子类型为电子,电子源为固定点源并设置点的坐标值,设置电子源的能量为单能,电子源单方向出射;
计数建模,计数粒子类型为光子,计数栅元号为用于计数的几何模型栅元号,统计光子通量,为当前计数添加辅助能量计数卡设置能群。
作为本发明的可选实施方式,本发明的一种X射线图像校正方法中,所述按照能谱划分步长将X射线能谱分布划分成若干个单能的子能谱中的能谱划分步长采用等步长划分方法确定;
可选地,所述的能谱步长为0.1keV-10keV范围。
本发明同时提供一种实现所述X射线图像校正方法的系统,包括:
几何及物理建模模块,建立不同贯穿长度L被测物体的X射线照射的几何模型、材料模型和源模型;
蒙卡模拟模块,基于物理建模开展蒙特卡罗粒子输运并统计计数计算投影值P;
多项式拟合模块,将不同贯穿长度L和投影值P进行多项式拟合;
图像校正模块,利用得到贯穿长度L和投影值P的拟合曲线校正实际的多能投影值Pm,得到理想的单能投影值Ps,从而应用到X射线图像的反投影重建。
本发明还提供一种计算机可读存储介质,存储有计算机可执行程序,所述计算机可执行程序被执行时,实现所述的一种X射线图像校正方法。
与现有技术相比,本发明的有益效果:
本发明的一种X射线图像校正方法,首先模拟X射线球管的X射线能谱分布,然后确定能谱划分步长将能谱划分为若干个单能X射束,利用蒙卡方法模拟划分步长后的单能X射束穿过不同厚度L被测物体到达探测器的光子通量,再把各个单能射束计算的光子通量转换成探测器平板响应并作加权、归一化处理,得到多能投影值Pm,最后利用多项式拟合二者的函数关系,取靠近原点的几组投影值进行以贯穿长度L为自变量的直线拟合,该直线的斜率k为等效衰减系数,由等效单能投影Ps滤波反投影重建图像。因此,本发明的一种X射线图像校正方法,基于蒙卡模拟克服了硬化校正对实验条件所限制的问题,具有灵活性高、模拟精度高的优势,能有效提高重建图像的均匀性。
附图说明:
图1本发明背景技术中单能和多能时,投影值与贯穿长度的函数曲线对比图;
图2本发明实施例一种X射线图像校正方法的流程图;
图3本发明实施例X射线能谱分布模拟几何建模;
图4本发明实施例光子通量分布模拟几何建模;
图5本发明实施例100kV管电压下投影值和贯穿长度多项式曲线;
图6本发明实施例120kV管电压下投影值和贯穿长度多项式曲线;
图7本发明实施例100kV条件下硬化校正前后均匀性检测图域中第250行灰度值曲线;
图8本发明实施例120kV条件下硬化校正前后均匀性检测图域中第250行灰度值曲线;
图9本发明实施例100kV条件下重建CBCT的图像均匀性横截面图;
图10本发明实施例120kV条件下重建CBCT的图像均匀性横截面图;
图11本发明实施例在CBCT图像上选取感兴趣区域的分布示例图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合附图,对本发明实施例中的技术方案进行清楚、完整的描述。显然,所描述的实施例是本发明的一部分实施例,而不是全部的实施例。
因此,以下对本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围,而是仅仅表示本发明的部分实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
需要说明的是,在不冲突的情况下,本发明中的实施例及实施例中的特征和技术方案可以相互组合。
应注意到:相似的标号和字母在下面的附图中表示类似项,因此,一旦某一项在一个附图中被定义,则在随后的附图中不需要对其进行进一步定义和解释。
在本发明的描述中,需要说明的是,术语“上”、“下”等指示的方位或位置关系为基于附图所示的方位或位置关系,或者是该发明产品使用时惯常摆放的方位或位置关系,或者是本领域技术人员惯常理解的方位或位置关系,这类术语仅是为了便于描述本发明和简化描述,而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作,因此不能理解为对本发明的限制。此外,术语“第一”、“第二”等仅用于区分描述,而不能理解为指示或暗示相对重要性。
参见图2所示,本实施例的一种X射线图像校正方法,包括:
模拟X射线球管管电压能量为E的X射线能谱分布;
按照能谱划分步长将X射线能谱分布划分成若干个单能的子能谱;
模拟各个子能谱X射线穿过不同贯穿长度L的被检测物体的光子通量;
将子能谱的光子通量转换成投影值P;
将贯穿长度L与各子能谱加权后的投影值P进行多项式拟合;
利用得到贯穿长度L和投影值P的拟合曲线校正实际的多能投影值Pm,得到理想的单能投影值Ps,实现X射线图像的硬化校正。
本实施例的一种X射线图像校正方法,首先模拟X射线球管的X射线能谱分布,然后确定能谱划分步长将能谱划分为若干个单能X射束,利用蒙卡方法模拟划分步长后的单能X射束穿过不同厚度L被测物体到达探测器的光子通量,再把各个单能射束计算的光子通量转换成探测器平板响应并作加权、归一化处理,得到多能投影值Pm,最后利用多项式拟合二者的函数关系,取靠近原点的几组投影值进行以贯穿长度L为自变量的直线拟合,该直线的斜率k为等效衰减系数,由等效单能投影Ps滤波反投影重建图像。因此,本实施例的一种X射线图像校正方法,基于蒙卡模拟克服了硬化校正对实验条件所限制的问题,具有灵活性高、模拟精度高的优势,能有效提高重建图像的均匀性。
作为本实施例的可选实施方式,本实施例的一种X射线图像校正方法中,所述模拟各个子能谱X射线穿过不同贯穿长度L的被检测物体的光子通量利用蒙卡模拟,包括:
光子通量分布模拟几何建模,创建表示光子源的第一几何体,创建表示被检测物体的第二几何体,创建表示探测器平板的第三几何体和第四几何体,所述的第二几何体、第三几何体和第四几何体依次排布在第一几何体的一侧,所述第四几何体布置在距离第一几何体Smax位置处,其中,所述Smax表征距离光子源无穷远位置,所述第二几何体的厚度设置表征被检测物体的贯穿长度L;
光子通量分布模拟物理建模,为第二几何体、第三几何体和第四几何体模型选择相应的材料,源建模设置源粒子类型为光子源,光子源为固定点源并设置点的坐标值,设置光子源的能量为单能,大小为X射线能谱分布划分成若干个单能的子能谱的大小,设置光子源的方向参数为各向异性,光子源单方向出射;
光子通量分布模拟计数建模,计数粒子类型为光子,计数栅元号为用于计数的几何模型栅元号,统计光子源单能射束穿过物体光子数。
具体地,参见图4所示,本实施例的光子通量分布模拟几何建模,选择Sphere工具创建一个r=2cm的小球(即第一几何体)表示光子源500,创建三个圆柱体(即第二几何体、第三几何体和第四几何体),一个表示为水层600,其余两个一个表示碳板700,一个表示表示碘化铯探测器平板400,最外侧的碘化铯探测器平板400放置在与光子源500距离Smax远处,铯探测器平板400与光子源500距离应尽量远,相对于其他几何体的距离成10倍以上,视为无穷远,目的在于忽略非正入射方向光子经散射照射到CsI感光材料的影响,图4中Smax示例取1000cm。本实施例水层600的厚度范围在0~20cm。
光子通量分布模拟几何建模用以建立光子感光场景,模拟光子源穿过水层、探测器平板表面、照射到CsI感光材料的物理路径,其中,光子源建模用于模拟光子源头,水层建模用于模拟不同厚度的被照射物体,探测器平板表面建模用于表示探测器平板的外表面,CsI建模用于表示光子的感光材料。
物理建模中为水层、探测器平板几何模型选择相应的材料,源建模设置源粒子类型为光子,源为固定点源并设置点的坐标值;设置源粒子的能量为单能,大小为计算出X射线能谱分布划分成若干个单能的子能谱的大小;设置源的方向参数为各向异性,参考方向(0,0,-1),源粒子单方向出射,出射方向与参考方向夹角的余弦值为1。
计数建模中计数粒子类型为光子,计数栅元号为用于计数的几何模型栅元号,统计子能谱穿过物体光子计数,选择t6类型计数,t6类型是蒙卡模拟软件的计数类型,表示统计光子的能量沉积。
作为本实施例的可选实施方式,本实施例的一种X射线图像校正方法,将利用蒙卡模拟的各个子能谱X射线穿过不同贯穿长度L的被检测物体的光子通量分别保存为单能通量,构建含有各单能通量的通量数据库,在进行X射线图像的硬化校正时,根据划分成的若干个单能的子能谱调用通量数据库中相应的单能通量,快速计算拟合曲线。
这样,本实施例的一种X射线图像校正方法,提出了单能的子能谱划分的方法,通过蒙卡模拟单能通量,保存为通量数据库备用,硬化校正时通过“能谱分布+子能谱划分+数据库”的方法,快速计算拟合曲线,应用于X射线图像硬化校正。对于给定能谱能快速计算拟合曲线,具有计算速度快、适应性高的优势,相比之下,现有技术对不同能量X射线都要进行蒙卡模拟得到拟合曲线,过程费时。
举例来说,若现有技术蒙卡模拟了170kV的拟合曲线,若要知道120kV的拟合曲线,还需要蒙卡模拟一次。本方法通过划分的单能的蒙卡模拟结果,直接推出不同能量下X射线的拟合曲线,省去了耗时的蒙卡模拟。
本实施例的一种X射线图像校正方法中,所述将子能谱的光子通量转换成投影值P包括:
将蒙卡模拟各单能X射线穿过被检测物体后得到光子通量,将单能射束模拟的光子通量转成探测器平板响应,利用探测器平板响应正比于射线强度,累加归一化后计算得到投影值,探测器平板响应计算公式为
其中,i表示能群,t表示厚度,μi表示第i个能群对应被检测物体的质量衰减系数,μen,i表示第i个能群对应CsI的质量能量吸收系数,Φi,0表示蒙卡模拟的第i个能群透过厚度为t的光子通量,表示第i个能群的平均能量,Ii,t表示第i个能群透过厚度为t的水层在CsI探测器吸收的能量,Ii,t正比于探测器平板响应,将i个能群的Ii,t进行累加做归一化处理得到It,表示加权后的多能X射束透过厚度为t的被检测物体在CsI测探测器吸收的能量,计算公式如下:
利用平板亮度计算投影值,公式如下:
本实施例的一种X射线图像校正方法中,所述将贯穿长度L与各子能谱加权后的投影值P进行多项式拟合包括:
采用投影值P作为自变量,材料的贯穿长度L作为因变量,进行N次多项式最小二乘法拟合,得到L与P的函数关系如下:
其中,N不小于4。
本实施例的一种X射线图像校正方法中,所述利用得到贯穿长度L和投影值P的拟合曲线校正实际的多能投影值Pm,得到理想的单能投影值Ps,采用如下公式计算:
Ps=kf-1(Pm)
作为CT或CBCT等图像的多项式硬化校正输入。
作为本实施例的可选实施方式,本实施例的一种X射线图像校正方法中,所述模拟X射线球管管电压能量为E的X射线能谱分布利用蒙卡模拟,包括:
几何建模,创建表示电子源的第五几何体,创建表示钨靶的第六几何体,创建表示铝滤过板的第七几何体,创建表示探测器平板的第八几何体,所述的第六几何体具有靶角θ的倾斜靶面,所述的第五几何体位于第六几何体的倾斜靶面的竖直上侧,所述的第七几何体、第八几何体依次布置在第六几何体的倾斜靶面的水平右侧;
物理建模,针对本模拟实验场景使用的第六几何体、第七几何体、第八几何体选择相应的材料进行材料建模,源建模设置源粒子类型为电子,电子源为固定点源并设置点的坐标值,设置电子源的能量为单能,电子源单方向出射;
计数建模,计数粒子类型为光子,计数栅元号为用于计数的几何模型栅元号,选择t4类型计数,t4类型是蒙卡模拟软件的计数类型,t4表示统计光子通量,为当前计数添加辅助能量计数卡设置能群。
具体地,参见图3所示,本实施例几何建模首先选择Sphere工具创建一个r=2cm的小球(即第五几何体)表示电子源100,创建一个圆柱为底、圆锥为顶的钨靶200(即第六几何体),其中圆锥的高h、底面半径为r,tanθ=h/r,θ为靶角,最后创建两个长方体(即第七几何体和第八几何体),一个表示为铝滤过板300,一个表示碘化铯探测器平板400。本实施例铝滤过板300的厚度为4.5mm,碘化铯探测器平板400与钨靶200之间的距离S1为750mm。
物理建模首先进行材料建模,对本模拟实验场景使用的钨靶、铝滤过、碘化铯探测器选择相应的材料。源建模中设置源粒子类型为电子,源为固定点源并设置点的坐标值;设置源粒子的能量为单能;设置源的方向参数为各向异性,参考方向(0,0,-1),源粒子单方向出射,出射方向与参考方向夹角的余弦值为1。
计数建模中计数粒子类型为光子,计数栅元号为用于计数的几何模型栅元号,选择t4类型计数,为当前计数添加辅助能量计数卡设置能群。
本实施例的一种X射线图像校正方法中,所述按照能谱划分步长将X射线能谱分布划分成若干个单能的子能谱中的能谱划分步长采用等步长划分方法确定。
可选地,所述的能谱步长为0.1keV-10keV范围。
本实施例同时提供一种实现所述X射线图像校正方法的系统,包括:
几何及物理建模模块,建立不同贯穿长度L被测物体的X射线照射的几何模型、材料模型和源模型;
蒙卡模拟模块,基于物理建模开展蒙特卡罗粒子输运并统计计数计算投影值P;
多项式拟合模块,将不同贯穿长度L和投影值P进行多项式拟合;
图像校正模块,利用得到贯穿长度L和投影值P的拟合曲线校正实际的多能投影值Pm,得到理想的单能投影值Ps,从而应用到X射线图像的反投影重建。
为验证本发明方法,利用最小二乘法进行多能投影值P为自变量,贯穿长度L为因变量的四阶多项式拟合,得到L=f(P)曲线,分别研究了100kV、120kV X射线能谱划分步长为0.5keV的情况,得到的拟合曲线如图5和图6所示。并应用到投影图像重建CBCT,对重建后的CBCT进行了影像均匀性的评价。均匀性评价CBCT图像在成像区域内重现特定灰度值的能力,越小越好。
在均匀性测试对应模块,取横断面图像区域中第250行灰度值进行均匀性比较,绘制灰度值相对于像素位置关系曲线,如图7、图8所示。由二者曲线可知,由于多能谱X射线硬化效应,使得图像出现杯状伪影,即贯穿长度大的区域暗,且随着贯穿长度增大亮度越来越暗,圆柱中心区域的亮度明显低于圆柱边缘的亮度。经过硬化校正后,杯状伪影得到了明显的改善,圆柱亮度值更加均匀,反映在曲线上为校正后曲线更加平坦。
在图像区域选取5个不同的位置画直径为1cm的圆,得到5个区域,5个区域分布如图11所示,5个区域圆的直径都为1cm,其中区域E为中心感兴趣区域,区域A-D为外部均匀分布地的感兴趣区域。根据CBCT图像均匀性的标准要求,中心感兴趣区域的参考HU值与外部4个感兴趣区域平均HU值之差的绝对值的最大值≤30HU。
如图9和图10,对比了100kV和120kV条件下重建CBCT的均匀性,并与CBCT图像均匀性的标准进行对比。从表可以看出,经硬化校正后重建的CBCT影像均匀性优于未校正的CBCT均匀性,且达到了均匀性的标准要求,硬化校正后提高了CBCT图像质量。
表1 100kV条件下硬化校正前后重建CBCT的均匀性比较
表2 120kV条件下硬化校正前后重建CBCT的均匀性比较
本实施例还提供一种计算机可读存储介质,存储有计算机可执行程序,所述计算机可执行程序被执行时,实现所述的一种X射线图像校正方法。
本实施例所述计算机可读存储介质可以包括在基带中或者作为载波一部分传播的数据信号,其中承载了可读程序代码。这种传播的数据信号可以采用多种形式,包括但不限于电磁信号、光信号或上述的任意合适的组合。计算机可读存储介质还可以是可读存储介质以外的任何可读介质,该计算机可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。计算机可读存储介质上包含的程序代码可以用任何适当的介质传输,包括但不限于无线、有线、光缆、RF等等,或者上述的任意合适的组合。
本实施例还提供了一种电子设备,包括处理器和存储器,所述存储器用于存储计算机可执行程序,当所述计算机程序被所述处理器执行时,所述处理器执行所述的一种X射线图像校正方法。
电子设备以通用计算设备的形式表现。其中处理器可以是一个,也可以是多个并且协同工作。本发明也不排除进行分布式处理,即处理器可以分散在不同的实体设备中。本发明的电子设备并不限于单一实体,也可以是多个实体设备的总和。
所述存储器存储有计算机可执行程序,通常是机器可读的代码。所述计算机可读程序可以被所述处理器执行,以使得电子设备能够执行本发明的方法,或者方法中的至少部分步骤。
所述存储器包括易失性存储器,例如随机存取存储单元(RAM)和/或高速缓存存储单元,还可以是非易失性存储器,如只读存储单元(ROM)。
应当理解,本发明的电子设备中还可以包括上述示例中未示出的元件或组件。例如,有些电子设备中还包括有显示屏等显示单元,有些电子设备还包括人机交互元件,例如按钮、键盘等。只要该电子设备能够执行存储器中的计算机可读程序以实现本发明方法或方法的至少部分步骤,均可认为是本发明所涵盖的电子设备。
通过以上对实施方式的描述,本领域的技术人员易于理解,本发明可以由能够执行特定计算机程序的硬件来实现,例如本发明的系统,以及系统中包含的电子处理单元、服务器、客户端、手机、控制单元、处理器等。本发明也可以由执行本发明的方法的计算机软件来实现,例如由微处理器、电子控制单元,客户端、服务器端等执行的控制软件来实现。但需要说明的是,执行本发明的方法的计算机软件并不限于由一个或特定个的硬件实体中执行,其也可以是由不特定具体硬件的以分布式的方式来实现。对于计算机软件,软件产品可以存储在一个计算机可读的存储介质(可以是CD-ROM,U盘,移动硬盘等)中,也可以分布式存储于网络上,只要其能使得电子设备执行根据本发明的方法。
以上实施例仅用以说明本发明而并非限制本发明所描述的技术方案,尽管本说明书参照上述的各个实施例对本发明已进行了详细的说明,但本发明不局限于上述具体实施方式,因此任何对本发明进行修改或等同替换;而一切不脱离发明的精神和范围的技术方案及其改进,其均涵盖在本发明的权利要求范围当中。
Claims (10)
1.一种X射线图像校正方法,其特征在于,包括:
模拟X射线球管管电压能量为E的X射线能谱分布;
按照能谱划分步长将X射线能谱分布划分成若干个单能的子能谱;
模拟各个子能谱X射线穿过不同贯穿长度L的被检测物体的光子通量;
将子能谱的光子通量转换成投影值P;
将贯穿长度L与各子能谱加权后的投影值P进行多项式拟合;
利用得到贯穿长度L和投影值P的拟合曲线校正实际的多能投影值Pm,得到理想的单能投影值Ps,实现X射线图像的硬化校正。
2.根据权利要求1所述的一种X射线图像校正方法,其特征在于,所述模拟各个子能谱X射线穿过不同贯穿长度L的被检测物体的光子通量利用蒙卡模拟,包括:
光子通量分布模拟几何建模,创建表示光子源的第一几何体,创建表示被检测物体的第二几何体,创建表示探测器平板的第三几何体和第四几何体,所述的第二几何体、第三几何体和第四几何体依次排布在第一几何体的一侧,所述第四几何体布置在距离第一几何体Smax位置处,其中,所述Smax表征距离光子源无穷远位置,所述第二几何体的厚度设置表征被检测物体的贯穿长度L;
光子通量分布模拟物理建模,为第二几何体、第三几何体和第四几何体模型选择相应的材料,源建模设置源粒子类型为光子源,光子源为固定点源并设置点的坐标值,设置光子源的能量为单能,大小为X射线能谱分布划分成若干个单能的子能谱的大小,设置光子源的方向参数为各向异性,光子源单方向出射;
光子通量分布模拟计数建模,计数粒子类型为光子,计数栅元号为用于计数的几何模型栅元号,统计光子源单能射束穿过物体光子数。
3.根据权利要求2所述的一种X射线图像校正方法,其特征在于,将利用蒙卡模拟的各个子能谱X射线穿过不同贯穿长度L的被检测物体的光子通量分别保存为单能通量,构建含有各单能通量的通量数据库,在进行X射线图像的硬化校正时,根据划分成的若干个单能的子能谱调用通量数据库中相应的单能通量,快速计算拟合曲线。
4.根据权利要求2所述的一种X射线图像校正方法,其特征在于,所述将子能谱的光子通量转换成投影值P包括:
将蒙卡模拟各单能X射线穿过被检测物体后得到光子通量,将单能射束模拟的光子通量转成探测器平板响应,利用探测器平板响应正比于射线强度,累加归一化后计算得到投影值,探测器平板响应计算公式为
其中,i表示能群,t表示厚度,μi表示第i个能群对应被检测物体的质量衰减系数,μen,i表示第i个能群对应CsI的质量能量吸收系数,Φi,0表示蒙卡模拟的第i个能群透过厚度为t的光子通量,表示第i个能群的平均能量,Ii,t表示第i个能群透过厚度为t的水层在CsI探测器吸收的能量,Ii,t正比于探测器平板响应,将i个能群的Ii,t进行累加做归一化处理得到It,表示加权后的多能X射束透过厚度为t的被检测物体在CsI测探测器吸收的能量,计算公式如下:
利用平板亮度计算投影值,公式如下:
5.根据权利要求4所述的一种X射线图像校正方法,其特征在于,所述将贯穿长度L与各子能谱加权后的投影值P进行多项式拟合包括:
采用投影值P作为自变量,材料的贯穿长度L作为因变量,进行N次多项式最小二乘法拟合,得到L与P的函数关系如下:
其中,N不小于4。
6.根据权利要求1所述的一种X射线图像校正方法,其特征在于,所述利用得到贯穿长度L和投影值P的拟合曲线校正实际的多能投影值Pm,得到理想的单能投影值Ps,采用如下公式计算:
Ps=kf-1(Pm)
作为CT或CBCT等图像的多项式硬化校正输入。
7.根据权利要求1所述的一种X射线图像校正方法,其特征在于,所述模拟X射线球管管电压能量为E的X射线能谱分布利用蒙卡模拟,包括:
几何建模,创建表示电子源的第五几何体,创建表示钨靶的第六几何体,创建表示铝滤过板的第七几何体,创建表示探测器平板的第八几何体,所述的第六几何体具有靶角θ的倾斜靶面,所述的第五几何体位于第六几何体的倾斜靶面的竖直上侧,所述的第七几何体、第八几何体依次布置在第六几何体的倾斜靶面的水平右侧;
物理建模,针对本模拟实验场景使用的第六几何体、第七几何体、第八几何体选择相应的材料进行材料建模,源建模设置源粒子类型为电子,电子源为固定点源并设置点的坐标值,设置电子源的能量为单能,电子源单方向出射;
计数建模,计数粒子类型为光子,计数栅元号为用于计数的几何模型栅元号,统计光子通量,为当前计数添加辅助能量计数卡设置能群。
8.根据权利要求1所述的一种X射线图像校正方法,其特征在于,所述按照能谱划分步长将X射线能谱分布划分成若干个单能的子能谱中的能谱划分步长采用等步长划分方法确定;
可选地,所述的能谱步长为0.1keV-10keV范围。
9.一种实现如权利要求1-8任意一项所述X射线图像校正方法的系统,其特征在于,包括:
几何及物理建模模块,建立不同贯穿长度L被测物体的X射线照射的几何模型、材料模型和源模型;
蒙卡模拟模块,基于物理建模开展蒙特卡罗粒子输运并统计计数计算投影值P;
多项式拟合模块,将不同贯穿长度L和投影值P进行多项式拟合;
图像校正模块,利用得到贯穿长度L和投影值P的拟合曲线校正实际的多能投影值Pm,得到理想的单能投影值Ps,从而应用到X射线图像的反投影重建。
10.一种计算机可读存储介质,其特征在于,存储有计算机可执行程序,所述计算机可执行程序被执行时,实现如权利要求1-8任意一项所述的一种X射线图像校正方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311682374.1A CN117665900B (zh) | 2023-12-07 | 2023-12-07 | 一种x射线图像校正方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311682374.1A CN117665900B (zh) | 2023-12-07 | 2023-12-07 | 一种x射线图像校正方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117665900A true CN117665900A (zh) | 2024-03-08 |
| CN117665900B CN117665900B (zh) | 2025-09-02 |
Family
ID=90073047
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311682374.1A Active CN117665900B (zh) | 2023-12-07 | 2023-12-07 | 一种x射线图像校正方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117665900B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN119151837A (zh) * | 2024-08-21 | 2024-12-17 | 北京富通康影科技有限公司 | 一种快速光子计数ct投影数据校正方法 |
| CN119693485A (zh) * | 2024-12-04 | 2025-03-25 | 东南大学 | 一种ct图像环形伪影数据模拟及移除方法 |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102609908A (zh) * | 2012-01-13 | 2012-07-25 | 中国人民解放军信息工程大学 | 基于基图像tv模型的ct射束硬化校正方法 |
| US20140218362A1 (en) * | 2013-02-05 | 2014-08-07 | Carestream Health, Inc. | Monte carlo modeling of field angle-dependent spectra for radiographic imaging systems |
| US20180235562A1 (en) * | 2017-02-17 | 2018-08-23 | Toshiba Medical Systems Corporation | Combined sinogram- and image-domain material decomposition for spectral computed tomography (ct) |
| CN109919868A (zh) * | 2019-02-27 | 2019-06-21 | 西北工业大学 | 一种锥束ct射束硬化曲线侦测及投影加权校正方法 |
| CN110811660A (zh) * | 2019-10-25 | 2020-02-21 | 赛诺威盛科技(北京)有限公司 | 一种校正ct射线束硬化伪影的方法 |
| CN111487265A (zh) * | 2020-05-25 | 2020-08-04 | 中国人民解放军战略支援部队信息工程大学 | 一种结合投影一致性的锥束ct硬化伪影校正方法 |
| CN111707688A (zh) * | 2020-06-18 | 2020-09-25 | 清华大学 | 光子计数能谱ct成像中自适应能谱优化方法及其应用 |
| CN115078421A (zh) * | 2022-06-13 | 2022-09-20 | 中国工程物理研究院材料研究所 | 一种基于能谱分析的ct射线硬化伪影校正方法及系统 |
-
2023
- 2023-12-07 CN CN202311682374.1A patent/CN117665900B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102609908A (zh) * | 2012-01-13 | 2012-07-25 | 中国人民解放军信息工程大学 | 基于基图像tv模型的ct射束硬化校正方法 |
| US20140218362A1 (en) * | 2013-02-05 | 2014-08-07 | Carestream Health, Inc. | Monte carlo modeling of field angle-dependent spectra for radiographic imaging systems |
| US20180235562A1 (en) * | 2017-02-17 | 2018-08-23 | Toshiba Medical Systems Corporation | Combined sinogram- and image-domain material decomposition for spectral computed tomography (ct) |
| CN109919868A (zh) * | 2019-02-27 | 2019-06-21 | 西北工业大学 | 一种锥束ct射束硬化曲线侦测及投影加权校正方法 |
| CN110811660A (zh) * | 2019-10-25 | 2020-02-21 | 赛诺威盛科技(北京)有限公司 | 一种校正ct射线束硬化伪影的方法 |
| CN111487265A (zh) * | 2020-05-25 | 2020-08-04 | 中国人民解放军战略支援部队信息工程大学 | 一种结合投影一致性的锥束ct硬化伪影校正方法 |
| CN111707688A (zh) * | 2020-06-18 | 2020-09-25 | 清华大学 | 光子计数能谱ct成像中自适应能谱优化方法及其应用 |
| CN115078421A (zh) * | 2022-06-13 | 2022-09-20 | 中国工程物理研究院材料研究所 | 一种基于能谱分析的ct射线硬化伪影校正方法及系统 |
Non-Patent Citations (3)
| Title |
|---|
| 张益海 等: "锥束工业CT射束硬化校正方法", 无损检测, vol. 39, no. 6, 31 December 2017 (2017-12-31), pages 8 - 12 * |
| 曾钢 等: "基于蒙特卡罗模拟的射束硬化校正方法", 高能物理与核物理, vol. 30, no. 2, 15 February 2006 (2006-02-15), pages 178 - 182 * |
| 陈慧娟 等: "一种基于多能统计的射束硬化校正方法", CT理论与应用研究, vol. 19, no. 1, 31 March 2010 (2010-03-31), pages 21 - 27 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN119151837A (zh) * | 2024-08-21 | 2024-12-17 | 北京富通康影科技有限公司 | 一种快速光子计数ct投影数据校正方法 |
| CN119693485A (zh) * | 2024-12-04 | 2025-03-25 | 东南大学 | 一种ct图像环形伪影数据模拟及移除方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117665900B (zh) | 2025-09-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US10302578B2 (en) | Method for performing material decomposition using a dual-energy X-ray CT and corresponding dual-energy X-ray CT apparatus | |
| CN107802280B (zh) | 校正曲线生成方法、投影图像的校正方法、系统及存储介质 | |
| US9916656B2 (en) | Method for processing radiographic image and radiography system | |
| CN117665900A (zh) | 一种x射线图像校正方法及系统 | |
| CN101566590B (zh) | 面阵探测器射线数字成像中的散射强度分布获取方法 | |
| Dhaene et al. | A realistic projection simulator for laboratory based X-ray micro-CT | |
| US9589373B2 (en) | Monte carlo modeling of field angle-dependent spectra for radiographic imaging systems | |
| Lazos et al. | Impact of flat panel‐imager veiling glare on scatter‐estimation accuracy and image quality of a commercial on‐board cone‐beam CT imaging system | |
| CN107330952B (zh) | 电子计算机断层扫描图像数据处理系统 | |
| CN104161536A (zh) | 一种基于互补光栅的锥束ct散射校正方法及其装置 | |
| Miceli et al. | Monte Carlo simulations of a high-resolution X-ray CT system for industrial applications | |
| CN107202805B (zh) | 基于卷积核的锥束ct散射伪影校正方法 | |
| Shi et al. | GPU-accelerated Monte Carlo simulation of MV-CBCT | |
| US9615807B2 (en) | Systems and methods for improving image quality in cone beam computed tomography | |
| CN116563413B (zh) | 校正参数确定方法、图像重建方法、装置、设备及介质 | |
| US10448904B2 (en) | Decomposition method and apparatus based on basis material combination | |
| US10079078B2 (en) | Method for correcting a spectrum | |
| CN115078421A (zh) | 一种基于能谱分析的ct射线硬化伪影校正方法及系统 | |
| Danyk et al. | Using clustering analysis for determination of scattering kernels in X-ray imaging | |
| Yu et al. | Heel effect adaptive flat field correction of digital x‐ray detectors | |
| US11409019B1 (en) | Device for producing high resolution backscatter images | |
| JP7502868B2 (ja) | 画像処理装置および方法、プログラム | |
| Monnin et al. | Task-based detectability in anatomical background in digital mammography, digital breast tomosynthesis and synthetic mammography | |
| CN115330895A (zh) | 一种散射校正方法、装置、设备及存储介质 | |
| Doshi et al. | Characterization of an indirect X-ray imaging detector by simulation and experiment |
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 |