CN106725565B - 一种稀疏投影下的锥束xct成像质量评估方法 - Google Patents
一种稀疏投影下的锥束xct成像质量评估方法 Download PDFInfo
- Publication number
- CN106725565B CN106725565B CN201611031469.7A CN201611031469A CN106725565B CN 106725565 B CN106725565 B CN 106725565B CN 201611031469 A CN201611031469 A CN 201611031469A CN 106725565 B CN106725565 B CN 106725565B
- Authority
- CN
- China
- Prior art keywords
- image
- cone
- imaging
- noise
- noise power
- 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.)
- Expired - Fee Related
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 83
- 238000001303 quality assessment method Methods 0.000 title claims abstract description 18
- 238000000034 method Methods 0.000 claims abstract description 45
- 238000001228 spectrum Methods 0.000 claims abstract description 40
- 238000012546 transfer Methods 0.000 claims abstract description 17
- 238000005070 sampling Methods 0.000 claims abstract description 12
- 210000000746 body region Anatomy 0.000 claims abstract description 11
- 230000011218 segmentation Effects 0.000 claims abstract description 9
- 238000013441 quality evaluation Methods 0.000 claims abstract description 4
- 230000009466 transformation Effects 0.000 claims description 10
- 230000003595 spectral effect Effects 0.000 claims description 8
- 238000011156 evaluation Methods 0.000 claims description 7
- 238000010606 normalization Methods 0.000 claims description 4
- 230000008901 benefit Effects 0.000 claims description 3
- 239000000463 material Substances 0.000 claims description 3
- 238000009499 grossing Methods 0.000 claims description 2
- 229920000058 polyacrylate Polymers 0.000 claims description 2
- 238000011160 research Methods 0.000 description 6
- 238000013461 design Methods 0.000 description 4
- 239000011159 matrix material Substances 0.000 description 4
- 230000005855 radiation Effects 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 230000018109 developmental process Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 238000000691 measurement method Methods 0.000 description 3
- 238000010521 absorption reaction Methods 0.000 description 2
- 238000002059 diagnostic imaging Methods 0.000 description 2
- 238000011002 quantification Methods 0.000 description 2
- WFKWXMTUELFFGS-UHFFFAOYSA-N tungsten Chemical compound [W] WFKWXMTUELFFGS-UHFFFAOYSA-N 0.000 description 2
- 238000005311 autocorrelation function Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000003795 chemical substances by application Substances 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 238000013170 computed tomography imaging Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000001816 cooling Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 210000004872 soft tissue Anatomy 0.000 description 1
- 238000004611 spectroscopical analysis Methods 0.000 description 1
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/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
- A61B6/03—Computed tomography [CT]
- A61B6/032—Transmission computed tomography [CT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- 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/30—Subject of image; Context of image processing
- G06T2207/30168—Image quality inspection
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Medical Informatics (AREA)
- General Health & Medical Sciences (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Physics & Mathematics (AREA)
- Radiology & Medical Imaging (AREA)
- Theoretical Computer Science (AREA)
- Molecular Biology (AREA)
- Pulmonology (AREA)
- Pathology (AREA)
- High Energy & Nuclear Physics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Biophysics (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- Optics & Photonics (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Quality & Reliability (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Image Analysis (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明公开一种稀疏投影下的锥束XCT成像质量评估方法,该方法在得到重建出仿体的三维图像后,进行成像***质量评估,包括:选取多个三维重建结果切片,沿z轴方向叠加平均三维重建结果切片得到二维图像,并得到极坐标系下图像;进行阈值分割,识别出目标体区域,然后通过质心法确定目标体区域的中心。设定过采样像素大小,将极坐标系下图像中各个同心圆环范围内的像素点进行像素值平均,而后按照离开圆心的距离排序,由此得到过采样边缘扩散函数曲线ESF。求调制传递函数;计算锥束XCT成像***的二维噪声功率谱。得到一维噪声功率谱,计算噪声等效量子数NEQ。
Description
技术领域
本发明属于生物医学工程及计算机领域,涉及一种锥束XCT成像质量评估方法。
背景技术
近年来,全球经济发展速度极快,人均寿命也在不断提升,人们对于医疗的依赖程度也在不断地加大,同时人们对于各种无创医疗手段也更加追求。X线设备作为医疗影像学中重要的仪器设备之一,自然成为了人们所关注的焦点。在提升X线设备成像质量的同时,降低X线照射剂量则是广大患者和医疗研究人员追求的目标之一。如何用最小的X线照射剂量,产生质量清晰的影像,是X线设备研究人员关注的重点。目前,由于锥束XCT具有成像速度快,射线剂量相对较低等优点,在临床应用中发展迅速,而低剂量锥束XCT研究的关注也愈发广泛。降低锥束XCT射线剂量主要有两种方式:一是改变管电流和管电压,这样的结果是图像噪声较大,图像的信噪比低;二是减少投影次数,但这种情况下采用传统的重建方法会造成重建时数据不一致性提高,重建质量下降。
针对上述问题,在近年来的医学影像领域,发展稀疏投影下的锥束XCT成像方法以保证重建图像质量的研究得到了广泛关注,这种稀疏投影为前提的XCT重建方法,由于减少了投影成像数量,不仅减少辐射剂量,减少扫描时间,也避免了因为长时间扫描中病人的移动而造成的运动伪影,因此对于新一代锥束XCT***的发展具有重要的意义。但到目前为止,这类方法尚未得到临床的尝试和应用,其中一个重要的原因在于缺乏合理可靠的成像质量客观定量评估方法,目前相关的文献中只是给出稀疏投影锥束XCT重建与传统方法的定性比较,缺乏定量化标准。
对于稀疏投影锥束XCT图像重建而言,需要在三维成像空间内进行成像***和成像方法的评估,比传统的投影成像***,例如X射线透视这种二维成像***的性能评估更加复杂。到目前为止,尚无被广泛认可的锥束XCT成像***性能评估的标准,更没有针对稀疏投影下锥束XCT成像质量的评估。
当前国际上公认的成像***性能评价参数包括调制传递函数(ModulationTransfer Function,MTF),噪声功率谱(Noise Power of Spectra,NPS),和噪声等效量子数(Noise Equivalent Quanta,NEQ)等,但仅给出了针对X射线投影成像***的标准,且对应的标定仿体也仅适合于二维投影***。而在三维的锥束XCT成像***下,成像***及方法的差别造成二维投影***下采用的标定仿体和方法不再适用,更不适用于稀疏投影下的锥束XCT成像质量评估。目前有部分学者提出采用直径小于像素值大小的钨丝来实现锥束XCT的调制传递函数测量,但是这种方法中,由于钨丝直径小,成像结果质量差,且这种方法不易实现过采样,导致所得到的MTF结果缺失信息甚至失真。
因此,为了实现稀疏投影下的锥束XCT成像质量评估,设计新型仿体,发展出一种有效的评估方法,实现成像质量评估的标准化及定量化,已经成为稀疏投影下锥束XCT成像质量评估的研究重点,这对于成像***的发展应用具有十分重要的意义。
发明内容
本发明的主旨是提供一种稀疏投影下的锥束XCT成像质量评估方法,以此解决锥束XCT成像中稀疏投影三维重建的成像质量评估的基本问题:实现稀疏投影下的锥束XCT成像结果的调制传递函数,噪声功率谱及噪声等效量子数的准确测量,为评估欠定条件下稀疏投影锥束XCT成像方法的性能提供有力条件。技术方案如下:
一种稀疏投影下的锥束XCT成像质量评估方法,该方法将仿体置于锥束XCT成像***的点光源下方的成像腔,成像腔在水平平面上绕着成像腔中轴旋转,每旋转一定角度,进行一次曝光成像,利用TV正则化方法重建出仿体的三维图像;基于此三维图像进行成像***质量评估,分为两个部分,
第一部分,计算稀疏投影下的锥束XCT成像***的调制传递函数MTF
①选取多个三维重建结果切片,沿z轴方向叠加平均三维重建结果切片得到二维图像p(x,y),(x,y)为图像像素在直角坐标系下位置,并将其从直角坐标系映射到极坐标系,得到极坐标系下图像p(r,θ)。
②对图像p(x,y)采用阈值分割将灰度图像p(x,y)映射为二值图像pnorm(x,y),这里为图像中空气部分的像素平均值,为图像中成像质量评估仿体部分的像素平均值,将仿体和空气部分分割,识别出目标体区域,然后通过质心法确定目标体区域的中心。
③设定过采样像素大小Δr,由此形成半径分别为Δr,2Δr,3Δr,4Δr……的一系列同心圆环,将图像p(r,θ)中各个同心圆环范围内的像素点进行像素值平均,而后按照离开圆心的距离排序,由此得到过采样边缘扩散函数曲线ESF。
④对ESF求微分以后得到线扩散函数LSF,对LSF加hann窗后进行傅里叶变换,所得结果取模并进行零频率归一化后即为所求的调制传递函数MTF(u),这里的u为频率;
第二部分,求取稀疏投影下的锥束XCT成像***的噪声功率谱NPS:
⑤对于所选取的三维重建结果切片,分别在每一个切片的圆形区域内的同一区域截取方形ROI,计算该方形ROI的平均像素值,而后将原ROI图像减去该平均像素值,以此获得各个切片的噪声图像。
⑥将每一幅噪声图像都分成互不重叠的几个子块ROI,而后分别对所有子块ROI计算噪声功率谱,计算噪声功率谱时采用的傅里叶变换的点数与计算MTF时采用的点数一致,再将所得的噪声功率谱叠加平均,得到锥束XCT成像***的二维噪声功率谱。
⑦从二维噪声功率谱图中心沿横轴向右至终点,取横轴上下各7行功率谱值,而后将这14行功率谱值叠加平均,得到锥束XCT成像***的一维噪声功率谱NPSd(u)。
⑧根据公式计算得到噪声等效量子数NEQ,即得到成像质量空间分辨率和噪声水平的综合评价结果。
本发明的稀疏投影下的锥束XCT成像质量评估方法,该发明针对目前稀疏投影下锥束XCT成像缺乏定量评估方法的问题,利用圆柱体仿体,并设计相应的极坐标系下的过采样调制传递函数测量方法以及噪声功率谱测量方法,进而获取噪声等效量子数,以实现稀疏投影三维重建的成像质量的准确定量评估。该发明为准确获取稀疏投影锥束XCT三维重建的成像***的噪声等效量子数,有效评估其成像性能奠定了方法学基础,从而为稀疏投影锥束XCT三维重建方法在影像学临床实践的应用推广提供有力支持。
附图说明
图1成像***评估仿体设计
图2边缘扩散函数(ESF)形成示意图
图3噪声功率谱(NPS)计算中ROI选取示意图
图4锥束XCT成像***示意图
图5(a)(b)(c)分别为.不同投影角度数情况下,滤波反投影算法和TV正则化算法的MTF评估结果对比
图6(a)(b)(c)分别为不同投影角度数情况下,滤波反投影算法和TV正则化算法的NPS评估结果对比
图7.(a)(b)(c)分别为不同投影角度数情况下,滤波反投影算法和TV正则化算法的NEQ评估结果对比
具体实施方式
针对目前在稀疏投影锥束XCT三维重建方法定量评估方面的问题,本发明提出一种稀疏投影下的锥束XCT成像质量评估方法,可以准确获得与成像质量相关的过采样调制传递函数,噪声功率谱以及噪声等效量子数。最后得到的技术方案如下:
1.设置稀疏投影下旋转步进角度(比如每隔2度旋转物体一次,则对应180个投影),将设计的成像质量评估仿体安置于成像腔中,而后在物体旋转到设定角度时,进行曝光成像;
2.利用TV正则化方法重建出仿体的三维图像;
3.按照如下步骤计算稀疏投影下成像***的调制传递函数MTF:
①.沿z轴方向叠加平均三维重建结果切片(选取6个切片叠加)得到二维图像p(x,y),(x,y)为图像像素在直角坐标系下位置,并将其从直角坐标系映射到极坐标系,得到极坐标系下图像p(r,θ)。
②.对图像p(x,y)采用阈值分割将图像像素值映射到[0,1],即将灰度图像p(x,y)映射为二值图像pnorm(x,y),这里为图像中空气部分的像素平均值,为图像中成像质量评估仿体(圆形目标体)部分的像素平均值,将物体和空气部分分割,识别出目标体区域,然后通过质心法确定圆形目标体区域的中心。
③.设定过采样像素大小Δr,由此形成半径分别为Δr,2Δr,3Δr,4Δr……的一系列同心圆环,将图像p(r,θ)中各个同心圆环范围内的像素点进行像素值平均,而后按照离开圆心的距离排序,由此即可得到过采样边缘扩散函数曲线(ESF)。
④.对ESF求微分以后得到LSF,对LSF加hann窗后进行傅里叶变换,所得结果取模并进行零频率归一化后即为所求的MTF(u),这里的u为频率。
4.按照如下步骤计算稀疏投影下成像***的噪声功率谱NPS:
①.选择计算MTF时采用的6个三维重建结果切片,而后分别在每一个切片的圆形区域内的同一区域截取方形ROI,计算该方形ROI的平均像素值,而后将原ROI图像减去该平均像素值,以此获得6幅噪声图像。
②.将每一幅噪声图像都分成互不重叠的4个方形ROI,而后分别对所有子块ROI计算噪声功率谱,计算噪声功率谱时采用的傅里叶变换的点数与计算MTF时采用的点数一致。再将所得的噪声功率谱叠加平均,可以得到锥束XCT成像***的二维噪声功率谱。
③.从二维噪声功率谱图中心沿横轴向右至终点,取横轴上下各7行功率谱值,而后将这14行功率谱值叠加平均,得到锥束XCT成像***的一维噪声功率谱NPSd(u)。
5.根据公式计算得到噪声等效量子数NEQ,即得到成像质量空间分辨率和噪声水平的综合评价结果。
下面对本发明进行详细说明。
1稀疏投影下的锥束XCT成像方法
对于锥束CT成像***,成像模型可以用下式来表示:
其中I是到达探测器的光子数,I0是入射的光子数,l为从X射线源到探测器的积分直线,即X射线的路径,μ为待求目标区域吸收系数分布。实际采样和重建过程均是离散化模型,第i条射线穿过物体之后,到达探测器的光子数为:
其中,投影矩阵M的元素mi,j表示第i条射线穿过图像第j个像素点的衰减系数。
得到的第i条射线的投影数据pi可以表示为:
其中pi为第i条射线对应的投影数据。那么将上式表示为矩阵形式,数字锥束XCT***可以用如下模型描述:
P=Mμ (4)
其中P为图像投影,M为CT***投影矩阵。
对于稀疏投影下的锥束XCT投影数据重建图像问题,从数学角度讲是求解欠定线性方程组的问题,即用投影数据P和***矩阵M求解原始图像μ。通过引入了待重建图像全变分(total variation,TV)先验信息,将上述问题等价于求解以下优化问题:
其中,第一项是保真项,为了满足成像模型,即公式(4)。第二项是基于对图像分段常数或者分段连续的假设,通过TV正则化方法引入分段常数的约束条件,锐化重建图像的边缘。||TV(·)||1为TV范数,实际上是图像梯度的L1-范数,τ为非负的正则化参数,用来平衡前述两项。另外,由于吸收系数μ是非负的,还需要在公式(3)满足的同时加入非负约束,即μ≥0。
2成像***评估仿体设计
为了满足稀疏投影下的锥束XCT成像质量的评估,设计一种适合用于成像***的调制传递函数(Modulation Transfer Function,MTF)及噪声功率谱(Noise PowerSpectrum)测量的仿体。仿体材料选择接近于人体软组织的聚丙烯酸酯,圆柱体边界平滑,内部材料均匀,所设计的仿体为直径20厘米,厚度4厘米的圆柱体,如图1所示。
通过改变投影数量,进行不同投影数量下的锥束XCT重建,以此实现基于稀疏投影的锥束XCT成像质量的准确评估。
3成像质量的调制传递函数(MTF)评估方法
利用所设计的仿体来测量MTF,采用的原理是求得径向过采样边缘扩散函数(EdgeSpread Function,ESF),然后求微分得到线扩散函数(Line Spread Function,LSF),对LSF进行傅里叶变换并取模求得MTF。基本的步骤如下:
1)沿z轴方向叠加平均三维重建结果切片(选取6个切片叠加)得到二维图像p(x,y),此时的二维图像可以认为是***的点扩散函数从三维降维到二维
这里(x,y)为图像像素在直角坐标系下位置,Δz是z轴叠加平均的切片数量。这里设定轴向分辨率是各向同性的,即沿x轴和y轴的空间分辨率相同。将二维图像p(x,y)从直角坐标系映射到极坐标系,得到极坐标系下图像p(r,θ),其中代表着轴向的空间变量,是对应的极坐标下的极角。
2)用阈值分割的方法识别目标体区域,然后通过质心法确定圆形目标体区域的中心:阈值分割时将图像值映射到[0,1]内,0代表空气,1代表仿体区域的平均值,如下式:
即将灰度图像p(x,y)映射为二值图像pnorm(x,y),这里为图像中空气部分的像素平均值,为图像中成像质量评估仿体(圆形目标体)部分的像素平均值。
3)设定过采样像素大小Δr(比如为原像素大小的1/10,则对应过采样倍数为10倍),由此形成半径分别为Δr,2Δr,3Δr,4Δr……的一系列同心圆环,如图2所示给出了半径从Δr~10Δr的10个同心圆环,用不同颜色表示。将图像p(r,θ)中各个同心圆环范围内的像素点进行像素值平均,而后按照离开圆心的距离排序,由此即可得到图2右下角对应的过采样边缘扩散函数曲线(ESF)。
4)对ESF求微分以后得到线扩散函数(LSF),对LSF加hann窗后进行傅里叶变换,所得结果取模并进行零频率归一化后即为所求的MTF(u),这里的u为频率。
4成像质量的噪声功率谱(NPS)评估方法
噪声功率谱被定义为空间频率域上随机信号在每一个频率位置的方差。根据维纳-辛钦理论,通过计算自相关函数的傅立叶变换就可以获得噪声功率谱。但在实际计算中,噪声功率谱通常通过对噪声密度图像的二维傅立叶变换取模来直接获取,具体计算公式如下:
这里FT[f(x,y)]是噪声图像f(x,y)的傅立叶变换,X和Y分别是图像在x方向和y方向的最大尺寸。
在数字成像***中,上述的噪声功率谱的连续表达形式通常需要转换为离散形式进行计算。因此,需要重新计算连续噪声图像函数f(x,y)在一系列离散位置点的像素值,这些位置点是x=mΔx,y=nΔy,m=0…M,n=0…N,这里X=MΔx,Y=NΔy。这样就可以把连续图像表示为离散形式f(m,n)。相应的,空间频率变量u和v的函数也需要在离散频率位置重新计算,这些位置包括u=kΔu,v=lΔv,k=0,±1,±2…±K,l=0,±1,±2…±L。最大空间频率即奈奎斯特采样频率在两个频率轴方向分别为U=1/2Δx及V=1/2Δy。
为了估计数字成像***地噪声功率谱,上述计算功率谱的连续表达式也需要转换为离散形式。根据上述定义,可以将公式(8)改写为:
这里FTd[f(m,n)]是数字图像f(m,n)的二维傅立叶变换。
仍然采用图1所示的仿体,实现噪声功率谱测量的基本步骤如下:
1.选择计算MTF时采用的6个三维重建结果切片,而后分别在每一个切片的圆形区域内的同一区域截取方形ROI(如图3所示),计算该方形ROI的平均像素值,而后将原ROI图像减去该平均像素值,以此获得6幅噪声图像。
2.将每一幅噪声图像都分成互不重叠的4个方形ROI,而后分别对所有子块ROI按照公式(9)计算噪声功率谱,计算噪声功率谱时采用的傅里叶变换的点数与计算MTF时采用的点数一致。再将所得的噪声功率谱叠加平均,可以得到锥束XCT成像***的二维噪声功率谱。
3.从二维噪声功率谱图中心沿横轴向右至终点,取横轴上下各7行功率谱值,而后将这14行功率谱值叠加平均,得到锥束XCT成像***的一维噪声功率谱NPSd(u)。
5成像质量的噪声等效量子数(NEQ)评估方法
对于二维图像,其二维NEQ可以计算为:
因为具有空间频率依赖性的NEQ参数用量子来表示(或更精确地,每单位面积的量子),该指标能够客观地比较在不同成像***获取的图像的质量,对于在实际医学影像仪器和工业应用中可以方便地比较配置如散射光栅、探测器探元、光源X射线束、剂量等的不同以及重建算法的不同而带来的图像的差异。
在一维情况下,NEQ可以表示为如下形式
本发明采用美国BIOPTICS公司生产的Pixarray 100小动物数字放射成像***构建一套锥束XCT成像***,采用如图4采用的成像方式,所设计仿体安置在圆柱形物体承载容器位置。
我们将分别进行了180个角度、90个角度以及60个角度下对成像质量评估仿体进行投影。将上述三种情况下,分别采用目前锥束XCT采用的传统的滤波反投影算法以及TV正则化方法进行三维重建,而后依据本发明的技术流程依次对重建结果的中心6个切片求解MTF、NPS以及NEQ。
对不同投影角度数情况下的不同算法MTF结果之间比较如图5所示,从图中可以发现,在三种不同投影角度数情况下,TV正则化方法的MTF值均高于滤波反投影算法,滤波反投影方法得到的MTF在空间频率为0.3mm-1时已经达到极限空间分辨率(10%MTF),而对于TV正则化方法,三种情况下的MTF曲线中,极限空间分辨率都不低于0.7mm-1,远高于传统的滤波反投影方法。
对不同投影角度数情况下的不同算法NPS结果之间比较如图6所示,从图中可以发现,在三种不同投影角度数情况下,TV正则化方法的NPS均远低于传统的滤波反投影算法。
图7显示在不同投影角度情况下,传统滤波反投影方法与TV正则化方法对应的***NEQ结果之间的比较。从图中可以明显看到,TV正则化方法获得了远高于传统滤波反投影方法的NEQ值,即表明了TV正则化方法在获取的图像的空间分辨率及噪声方面具有更优的性能。
最终结果表明,在进行数字X射线成像***的调制传递函数测量时,通过本发明的基于X射线成像***探测器特性的调制传递函数测量新方法,充分利用了探测器ESF曲线的单调性和凹凸性,可以获得较传统刀口测量方法更为准确的MTF曲线。该方法的应用,将为准确测量***调制传递函数,有效评估放射成像***性能,深入开展放射影像学临床实践和研究提供有力支持。
Claims (1)
1.一种稀疏投影下的锥束XCT成像质量评估方法,该方法将仿体置于锥束XCT成像***的点光源下方的成像腔,所述的仿体为由聚丙烯酸酯制成的边界平滑,内部材料均匀的圆柱体;成像腔在水平平面上绕着成像腔中轴旋转,每旋转一定角度,进行一次曝光成像,利用TV正则化方法重建出仿体的三维图像;基于此三维图像进行成像***质量评估,分为两个部分,其中,
第一部分:计算稀疏投影下的锥束XCT成像***的调制传递函数MTF
①选取多个三维重建结果切片,沿z轴方向叠加平均三维重建结果切片得到二维图像p(x,y),(x,y)为图像像素在直角坐标系下位置,并将其从直角坐标系映射到极坐标系,得到极坐标系下图像p(r,θ);
②对图像p(x,y)采用阈值分割将灰度图像p(x,y)映射为二值图像pnorm(x,y),这里为图像中空气部分的像素平均值,为图像中成像质量评估仿体部分的像素平均值,将仿体和空气部分分割,识别出目标体区域,然后通过质心法确定目标体区域的中心;
③设定过采样像素大小Δr,由此形成半径分别为Δr,2Δr,3Δr,4Δr……的一系列同心圆环,将图像p(r,θ)中各个同心圆环范围内的像素点进行像素值平均,而后按照离开圆心的距离排序,由此得到过采样边缘扩散函数曲线ESF;
④对ESF求微分以后得到线扩散函数LSF,对LSF加hann窗后进行傅里叶变换,所得结果取模并进行零频率归一化后即为所求的调制传递函数MTF(u),这里的u为频率;
第二部分:求取稀疏投影下的锥束XCT成像***的噪声功率谱NPS:
⑤对于所选取的三维重建结果切片,分别在每一个切片的圆形区域内的同一区域截取方形ROI,计算该方形ROI的平均像素值,而后将原ROI图像减去该平均像素值,以此获得各个切片的噪声图像;
⑥将每一幅噪声图像都分成互不重叠的几个子块ROI,而后分别对所有子块ROI计算噪声功率谱,计算噪声功率谱时采用的傅里叶变换的点数与计算MTF时采用的傅里叶变换的点数一致,再将所得的噪声功率谱叠加平均,得到锥束XCT成像***的二维噪声功率谱;
⑦从二维噪声功率谱图中心沿横轴向右至终点,取横轴上下各7行功率谱值,而后将这14行功率谱值叠加平均,得到锥束XCT成像***的一维噪声功率谱NPSd(u);
⑧根据公式计算得到噪声等效量子数NEQ,即得到成像质量空间分辨率和噪声水平的综合评价结果。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201611031469.7A CN106725565B (zh) | 2016-11-18 | 2016-11-18 | 一种稀疏投影下的锥束xct成像质量评估方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201611031469.7A CN106725565B (zh) | 2016-11-18 | 2016-11-18 | 一种稀疏投影下的锥束xct成像质量评估方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN106725565A CN106725565A (zh) | 2017-05-31 |
CN106725565B true CN106725565B (zh) | 2019-07-19 |
Family
ID=58971671
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201611031469.7A Expired - Fee Related CN106725565B (zh) | 2016-11-18 | 2016-11-18 | 一种稀疏投影下的锥束xct成像质量评估方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN106725565B (zh) |
Families Citing this family (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109589126B (zh) * | 2018-10-18 | 2022-08-19 | 天津大学 | 基于宽光束小步长扫描方式的x射线发光断层成像方法 |
CN109668912B (zh) * | 2018-12-05 | 2022-01-04 | 马亮 | 伴床式x光机质量检验机构 |
CN109975335B (zh) * | 2019-03-07 | 2024-01-12 | 北京航星机器制造有限公司 | 一种ct检测方法及装置 |
JP7317655B2 (ja) * | 2019-09-26 | 2023-07-31 | 富士フイルムヘルスケア株式会社 | 医用画像処理装置および医用画像処理方法 |
CN110619680A (zh) * | 2019-10-23 | 2019-12-27 | 浙江大学深圳研究院 | 一种基于图变分的三维断层相位显微镜重建方法 |
CN111127546B (zh) * | 2019-11-25 | 2023-04-28 | 南京航空航天大学 | 一种基于极坐标变换的圆形目标中心定位方法及*** |
CN113436221A (zh) * | 2021-05-31 | 2021-09-24 | 华东师范大学 | 一种利用几何形状分层的图像分割弱标注方法 |
CN113298702B (zh) * | 2021-06-23 | 2023-08-04 | 重庆科技学院 | 基于大尺寸图像像素点的重新排序和分割方法 |
CN113822961B (zh) * | 2021-09-22 | 2024-04-26 | 广州博冠信息科技有限公司 | 对3d模型进行2d渲染的方法、装置、设备及介质 |
Family Cites Families (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6521886B2 (en) * | 2000-12-29 | 2003-02-18 | Ge Medical Systems Global Technology Company, Llc | Method of monitoring changes in the detective quantum efficiency of an x-ray detector |
CN102012466A (zh) * | 2010-10-28 | 2011-04-13 | 天津大学 | 数字x射线成像***的噪声测量方法 |
CN102204828B (zh) * | 2011-05-13 | 2013-04-10 | 天津大学 | 数字x射线成像***的调制传递函数精确测量方法 |
CN104287761A (zh) * | 2013-07-19 | 2015-01-21 | 南京普爱射线影像设备有限公司 | C形臂x射线装置以及测量mtf值的方法 |
CN103528840B (zh) * | 2013-09-29 | 2016-03-30 | 天津大学 | 基于x射线成像***探测器特性的调制传递函数测量方法 |
KR20150089382A (ko) * | 2014-01-27 | 2015-08-05 | 연세대학교 원주산학협력단 | 디지털 유방 엑스선 촬영용 시스템을 위한 검출양자효율의 검출방법, 검출시스템 및 그 방법을 실행하는 프로그램이 기록된 기록매체 |
CN104083177B (zh) * | 2014-07-08 | 2016-07-06 | 天津大学 | 数字x射线成像***的分辨率性能评估方法 |
CN105424731B (zh) * | 2015-11-04 | 2018-03-13 | 中国人民解放军第四军医大学 | 一种锥束ct的分辨率性能测量装置及标定方法 |
-
2016
- 2016-11-18 CN CN201611031469.7A patent/CN106725565B/zh not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
CN106725565A (zh) | 2017-05-31 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN106725565B (zh) | 一种稀疏投影下的锥束xct成像质量评估方法 | |
US8615118B2 (en) | Techniques for tomographic image by background subtraction | |
Li et al. | Refraction corrected transmission ultrasound computed tomography for application in breast imaging | |
US7142633B2 (en) | Enhanced X-ray imaging system and method | |
CN110807737A (zh) | 迭代图像重建框架 | |
US9949709B2 (en) | Techniques for suppression of motion artifacts in medical imaging | |
US8075183B2 (en) | Method of evaluating the resolution of a volumetric imaging system and image phantom used during the resolution evaluation | |
US20160217596A1 (en) | Joint reconstruction of electron density images | |
JP2016152916A (ja) | X線コンピュータ断層撮像装置及び医用画像処理装置 | |
US10679346B2 (en) | Systems and methods for capturing deep learning training data from imaging systems | |
CN114067013A (zh) | 用于经由单应重采样变换进行重投影和反投影的***和方法 | |
Herbst et al. | Dynamic detector offsets for field of view extension in C‐arm computed tomography with application to weight‐bearing imaging | |
Makeev et al. | Investigation of statistical iterative reconstruction for dedicated breast CT | |
Xia et al. | Dedicated breast computed tomography: Volume image denoising via a partial‐diffusion equation based technique | |
CN110428370A (zh) | 一种利用偏心旋转提高锥形束spect成像分辨率的方法 | |
Muckley et al. | Image reconstruction for interrupted-beam x-ray CT on diagnostic clinical scanners | |
Zhong et al. | A dual‐view digital tomosynthesis imaging technique for improved chest imaging | |
CN113167913A (zh) | 针对常规成像的光子计数的能量加权 | |
Ruiter et al. | 3d ultrasound computer tomography for breast cancer diagnosis at kit: an overview | |
CN103202704B (zh) | 半扫描位置的确定方法 | |
EP4182891A1 (en) | Deep learning for sliding window phase retrieval | |
O’Connor et al. | Development of an ensemble of digital breast object models | |
Miao et al. | Implementation of FDK reconstruction algorithm in cone-beam CT based on the 3D Shepp-Logan Model | |
Kong et al. | Evaluation of an analytic reconstruction method as a platform for spectral cone-beam CT | |
US12073491B2 (en) | Energy weighting of photon counts for conventional imaging |
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 | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190719 Termination date: 20211118 |