CN105631876A - 一种基于全局二值化的ct图像分辨率自动测试方法 - Google Patents

Abstract

一种基于全局二值化的CT图像分辨率自动测试方法，包括采用圆盘标准试件进行测试，获取圆盘标准试件的CT图像；对CT图像进行自动阈值分割的全局二值化处理，并且在全局二值化处理后的CT图像中找圆盘标准试件的圆心以及圆盘轮廓线，并且计算圆盘的半径。根据中找出的圆盘标准时间的圆心以及圆盘轮廓线和半径，计算CT图像的密度分辨率和空间分辨率。该基于全局二值化的CT图像分辨率自动测试方法，在图像的中心和轮廓的确定计算中利用全局二值化后的CT图像，使得圆心、半径和轮廓的确定更加准确，进而使得分辨率的测试也更加准确。

Description

一种基于全局二值化的CT图像分辨率自动测试方法

技术领域

本发明涉及一种CT图像分辨率的自动测试方法。

背景技术

CT技术不仅是放射诊断医学领域的里程碑，也是现代工业无损检测和勘探领域的重要手段。工业CT技术不受被检测物体材料、形状、表面状况等因素影响，能够给出被检测物体二维、三维直观图像，成为工业设备或部件无损检测和质量评估的重要手段。如加速器CT(能量＞1MV)在航天航空、石油勘测、汽车装配业中的应用，微焦点CT(能量低于100kV)在精密仪器制造、生物学等领域的应用，以及锥束CT在三维尺寸测量、逆向重构、考古等领域的应用，医学CT的发展和应用在医学领域更是不可或缺。

CT系统质量保证及图像可靠性评定是CT研究领域的重要方向，CT系统性能指标测试方法的研究已在医学、工业CT领域中开展。CT系统的空间分辨率和密度分辨率是重要性能表征参数，也是CT检测质量保证的关键因素。空间分辨率，指从CT图像中能够分辨特定的最小几何细节的能力，定量的表示为能分辨两个细节特征的最小间距，医学临床上体现为对小病灶或结构的成像能力，工业CT上体现为对细节特征(气孔、裂纹)的辨别能力；密度分辨率是CT系统分辨给定面积映射到CT图像上射线衰减系数差别(对比度)的能力，定量的表示为给定面积上能够分辨的细节(给定面积)与基体材料的最小对比度。

目前，国内工业CT设备生产厂家有清华固鸿、重庆真测、四川中物院等多家单位，国内仅加速器CT就多达60余套，而低能CT、微焦点CT据统计超过百套(多为进口、拼装设备)应用在国家国防、石油、汽车、铁道等重要行业，医学CT更是多种多样。仅在工业CT领域，不同厂家设备性能指标测试方法不一，测试标准试件种类繁多，通常采用标准物质(试件)的直接测试方法来测试CT系统的空间、密度分辨能力。如空间分辨率测试中采用具有周期性结构的标准试件进行扫描成像，分析CT图像中按一定规律排布的周期性结构图像(通常有线对、圆孔、条形孔等)，以视觉上能分辨单位距离内的最多条纹或圆孔数目来测定极限空间分辨率，密度分辨率测试中采用密度差试件是在一均匀的基体材料的特定范围含有其它种类材料，这些材料的密度与基体材料略有差异，通过测试一定范围的平均灰度值的差异，来区分密度差。该方法简单、直观，但受到人为主观因素影响。随着工业CT技术的发展，微CT、锥束CT的逐步应用，标准试件的制作的成本、难度加大。在医学CT中对密度分辨率要求较高，通常采用系列分布的低对比度试件制作的体模来测试密度分辨能力。

在CT关键性能参数自动测量技术发展方面，目前法国科学家TarrafTorfeh等人在研究中尝试将性能测试模拟软件做入医用CT设备中，通过输入设备参数即可实时计算出该参数下的各项性能指标，该项技术对医学CT性能指标实时监测中发挥重要作用，而在工业CT中，目前国内外设备厂家多采用标准试件及视觉分辨来测试系统空间、密度分辨率。

发明内容

本发明所要解决的技术问题是提供一种二维CT图像分辨率自动测试方法，该方法可以用于测试CT图像的空间分辨率和密度分辨率。

本发明解决上述技术问题所采用的技术方案为：一种基于全局二值化的CT图像分辨率自动测试方法，包括采用圆盘标准试件进行测试，其特征在于：包括以下步骤：

1)获取圆盘标准试件的CT图像；

2)对CT图像进行自动阈值分割的全局二值化处理，并且在全局二值化处理后的CT图像中找圆盘标准试件的圆心以及圆盘轮廓线，并且计算圆盘的半径。

3)根据步骤2)中找出的圆盘标准试件的圆心以及圆盘轮廓线和半径，计算CT图像的密度分辨率和空间分辨率。

为了使圆心、半径和圆周的确定更加准确，上述步骤2)中的自动阈值分割的全局二值化处理包括将步骤1)中测得的图像的灰度值分为0到m级，灰度值为i的像素数为n_i，则总的像素数为各灰度出现的概率为P_i＝n_i/N，将灰度值分为c₀＝{0-k}和c₁＝{k+1-m}，则c₀出现的概率为c₁出现的概率为 ${\mathrm{\ω}}_1=\underset{i=k+1}{\overset{m}{\Σ}}{P}_i=1-\mathrm{\ω}\left(k\right),$ 可以推出c₀平均值为 $u_0=\underset{i=0}{\overset{k}{\Σ}}{\mathrm{iP}}_i/{\mathrm{\ω}}_0=\frac{u\left(k\right)}{\mathrm{\ω}\left(k\right)},$ c₁的平均值为 $u_1=\underset{i=k+1}{\overset{m}{\Σ}}{\mathrm{iP}}_i/{\mathrm{\ω}}_1=\frac{u-u\left(k\right)}{1-\mathrm{\ω}\left(k\right)},$ 其中 $u_k=\underset{i=0}{\overset{k}{\Σ}}{\mathrm{iP}}_i$ 是阈值为k时的灰度平均值， $u=\underset{i=0}{\overset{m}{\Σ}}{\mathrm{iP}}_i$ 是整体图像的灰度平均值，这样两极之间的方差就可以用下式表示：

σ²(k)＝ω₀(u-u₀)²+ω₁(u₁-u)²

从0到m之间改变k值，计算上式求其最大值便是所求的阈值，然后根据阈值将CT图像二值化，背景为0，圆盘为1；

采用从图像四周向图像中心设置搜索线，查找最外侧的灰度值为1的边界点，作为圆盘的轮廓线；对于获得的圆盘轮廓坐标进行圆周拟合得出圆心，并且计算得出半径。

优选地，上述步骤3)中测试空间分辨率的方法具体包括先获取边界响应函数(ERF)曲线，然后根据先获取边界响应函数(ERF)曲线获得点扩散函数(PSF)曲线，最后根据点扩散函数(PSF)曲线，获得调制传递函数(MTF)曲线。

优选地，获取边界响应函数(ERF)曲线的步骤，具体包括根据步骤2)中确定的圆盘半径R和圆心，选取圆环区域，圆环区域的外径为R+15个像素，内径为R-15个像素，该圆环区域覆盖圆盘标准试件的CT图像的边界，计算圆环区域内所有像素点到圆心的距离，并将距离相同的像素点归为一组，计算每组像素点的灰度值的而平均值，建立每组像素点到圆心的距离与像素灰度值的平均值之间的关系曲线，即边界响应函数(ERF)曲线，并对该边界响应函数(ERF)曲线采用分段插值法和最小二乘立方拟合，获得平滑边界响应函数(ERF)曲线。

优选地，获取点扩散函数(PSF)曲线的具体步骤包括，对获取的边界响应函数(ERF)曲线进行分段拟合，并对每一次拟合得出的多项式进行求导，计算中间点在每个导数解析式对应的值，得出距离与到数值的关系函数，对函数进行归一化处理，得出点扩散函数(PSF)曲线。

优选地，该获取调制传递函数(MTF)曲线为将点扩散函数进行离散傅里叶变换后的幅值在零频处归一化得到。

优选地，步骤3)中测试密度分辨率的方法具体包括，根据步骤2)中获取的圆盘标准试件的圆心和半径，选取图像中心的方形数据分析区域，并且将方形数据分析区域进行n×n分割成方块，n可以为任意自然数，计算每个方块的像素灰度平均值然后通过每个方块内的像素灰度平均值计算整个ROI区域的标准偏差，即得出该n数下该方形数据分析区域的平均标准偏差σ，选取不同的n，将不同n下的平均标准偏差表示为占该方形数据区域的像素的灰度平均值的百分比，建立不同的方块尺寸与该百分比之间的关系曲线即对比度鉴别函数(CDF)曲线。

优选地，该圆盘标准试件的直径选择为工业CT系统成像范围尺寸的70％，厚度不小于20毫米。

与现有技术相比，本发明的优点在于该基于全局二值化的CT图像分辨率自动测试方法，能够快速、自动得测试出CT图像的空间分辨率和密度分辨率，并且由于在图像的中心和轮廓的确定计算中利用全局二值化后的CT图像，使得圆心、半径和轮廓的确定更加准确，进而使得分辨率的测试也更加准确。

附图说明

图1为本发明实施例的空间分辨率的测试方法的流程图。

图2为本发明实施例的空间分辨率的测试方法所用的圆盘的图像示意图。

图3为本发明实施例的空间分辨率的测试方法中得出的MTF曲线的示意图。

图4为本发明实施例的密度分辨率的测试方法的流程图。

图5为本发明实施例的密度分辨率的测试方法中的ROI区域的示意图。

图6为本发明实施例的密度分辨率的测试方法中得出的CDF曲线的示意图。

具体实施方式

以下结合附图实施例对发明作进一步详细描述。

本发明的CT图像分辨率的自动测试方法，采用一种高精度加工、材料密度均一的圆盘试件作为标准试件。一般圆盘标准试件直径选择为工业CT系统成像范围尺寸的70％，厚度不小于20毫米。

调制传递函数(MTF)曲线反映了系统对周期性细节响应的特性，是表征工业CT系统空间分辨率的重要方法，由于实际中理想的点状物体成像很难实现，因此可以采用圆盘边界计算调制传递函数，即，先测试计算圆盘的边界响应函数(ERF)，一阶求导后获得点扩散函数(PSF)，对点扩散函数进行傅里叶变换得到系统的调制传递函数(MTF)。

所以该测试CT图像的空间分辨率的方法包括以下步骤：如图1的流程图所示。

1、CT图像的获取。

将标准试件置于CT系统的转台中心位置，调节CT设备至正常工作状态，设置扫描参数，获取圆盘标准试件的CT图像。

2、图像初步处理阶段，在输入CT图像后，对整幅CT图像进行自动阈值分割的全局二值化处理，并且在全局二值化处理后的CT图像中找出圆盘标准试件的圆心以及圆盘轮廓线并且计算圆盘的半径。该方法具体包括：

自动阈值分割的全局二值化处理为采用最小二乘法的阈值分割。具体的实现方法如下：设一幅图像的灰度值可以分为0到m级，灰度值为i的像素数为n_i，则总的像素数为：

$N=\underset{i=0}{\overset{m}{\Σ}}{n}_i---\left(1\right)$

那么各灰度值出现的概率为P_i＝n_i/N，然后将灰度值分为两组：c₀＝{0-k}和c₁＝{k+1-m}，则c₀出现的概率为：

${\mathrm{\ω}}_0=\underset{i-0}{\overset{k}{\Σ}}{P}_i=\mathrm{\ω}\left(k\right)---\left(2\right)$

c₁出现的概率为：

${\mathrm{\ω}}_1=\underset{i=k+1}{\overset{m}{\Σ}}{P}_i=1-\mathrm{\ω}\left(k\right)---\left(3\right)$

由以上两式可以推出c₀平均值为：

$u_0=\underset{i=0}{\overset{k}{\Σ}}{\mathrm{iP}}_i/{\mathrm{\ω}}_0=\frac{u\left(k\right)}{\mathrm{\ω}\left(k\right)}---\left(4\right)$

c₁的平均值为：

$u_1=\underset{i=k+1}{\overset{m}{\Σ}}{\mathrm{iP}}_i/{\mathrm{\ω}}_1=\frac{u-u\left(k\right)}{1-\mathrm{\ω}\left(k\right)}---\left(5\right)$

其中是阈值为k时的灰度平均值，是整体图像的灰度平均值。这样两极之间的方差就可以用下式表示：

σ²(k)＝ω₀(u-u₀)²+ω₁(u₁-u)²(6)

从0到m之间改变k值，计算上式求其最大值便是所求的阈值。

阈值分割以后，进行CT图像二值化，背景为0，圆盘为1。采用从图像四周向图像中心设置搜索线，查找最外侧的灰度值为1的边界点，作为圆盘外环轮廓线。因此，即可准确地确定圆盘轮廓线。

如图2所示，然后，进行360度搜索，对于获得的圆盘轮廓坐标进行圆周拟合，计算得出圆心1，可以采用最小二乘法圆心拟合法，从圆盘轮廓中计算出圆心位置，然后再计算出圆盘的半径，任意选取过一条圆心的直线(一般选择平行于X轴的直线)，可进一步确定原始CT图像中灰度值分布情况。

3、获取边界响应函数(ERF)曲线，程序根据上述步骤中确定的圆盘半径R和圆心，选取一圆环区域2，该圆环区域的外径为R+15个像素、内径为R-15个像素，该圆环区域2涵盖圆盘CT图像的边界，如图2所示，图中实线方框为图像边界，圆形的实线部分为CT圆盘图像的轮廓，而虚线为选取的圆环区域2，对圆环区域内的像素灰度值进行处理。

具体为：计算圆环区域2内所有像素点到圆心1的距离，并将距离相同的像素点归为一组，即圆环区域2内同半径的像素点为一组，计算每组像素点的灰度值的平均值，建立距离与像素灰度值的平均值的关系曲线，即ERF曲线。

然后对ERF曲线进行插值，采用分段插值法，每组像素点插入不少于20个值，然后进行最小二乘立方拟合，建立距离与拟合像素灰度值的关系曲线，即平滑ERF曲线。采用这种获得ERF曲线的方法较传统方法更好的克服噪声的影响。

4、获取点扩散函数(PSF)曲线，在此数据处理阶段，对生成的边界响应函数(ERF)曲线进行分段拟合，并对每一次拟合得出的多项式进行求导，计算中间点在每个导数解析式对应的值，得出距离与到数值的关系函数，对函数进行归一化处理，得出点扩散函数(PSF)曲线。

具体为：设ERF为一维数组，数组长度为n，E_n为每个的值。m是多项式阶数为3，a是多项式系数为1。

设i＝0～n-20，从En中取20长度数组，起点为i，设为X[i]，进行多项式拟合，拟合公式如下：

$y\[i\]=\underset{j=0}{\overset{m}{\Σ}}{a}_j{(X\[i\])}^j$

对每次求得的y数组进行拟合以后求导，公式如下：

$f\left(t\right)=\frac{d}{dt}F\left(t\right)$

取求导后f数组的中间值f(10)作为该次求导的值，累计形成新的数组M，即为PSF曲线数组。

5、获取调制传递函数(MTF)曲线，调制传递函数函数是从点扩散函数进行离散傅里叶变换后的幅值在零频处归一化得到。

具体为：设PSF为数组X，其元素数量为N，进行FFT变换，公式如下

$Y=\underset{n=0}{\overset{N-1}{\Σ}}{X}_n{e}^{-j2\mathrm{\π}n/N},$ 其中n＝0,1,2……N-1。

Y为变换结果，取模|Y|。

取|Y|中最大值，Y/|Y|即为零频处归一化.

如图3所示，该图中的实线即为求得的MTF曲线，其中的点是用传统的线对卡法求得的MTF曲线的示意图。

而密度分辨率的测试，利用的原理为多个尺寸相同低对比度区域灰度值的平均值为随机变量，且服从高斯分布，则多个与低对比度区域尺寸相同的背景区域的CT平均值也服从高斯分布，且两个分布的标准偏差相同。两个分布的区别在于它们的期望平均值。以两个分布的中点为阈值，用以从背景中分离出低对比度物体，那么当两个分布的平均值相离3σ(σ是分布的标准偏差)时，超过阈值的背景分布曲线下的面积达到0.27％，同理，低对比度物体分布低于阈值的部分也是0.27％，即要以99.73％的置信度从背景中区分出低对比度物体，对比度需要为3σ。密度分辨率的测试步骤的流程图如图4所示。

具体步骤为：

1、CT图像数据获取，将标准试件置于转台中心位置，调节工业CT设备至正常工作状态，设置扫描参数，获取标准试件的CT图像。

2、图像初步处理阶段，在输入CT图像后，对整幅CT图像进行自动阈值的全局二值化方法处理，并且在全局二值化方法处理后的CT图像中找出标准试件圆盘的中心以及圆盘轮廓线并且计算圆盘的半径。该步骤的具体方法可以参见上述对空间分辨率的测试中的步骤2。

3、在圆盘CT图像中心区域特定范围内选择一方形数据分析区域，即ROI区域，该方形数据区域对角线长度大约是圆盘直径的60％，并且中心位于圆盘的圆心，对方形区域按照尺寸大小进行n×n分割，n可以为任意自然数，如图5所示，该图中示出的n为5，即将该方形数据分析区域分割成5×5个方块，计算每个方块的像素灰度平均值，然后通过每个方块内的像素灰度平均值计算整个ROI区域的标准偏差，即得出该n数下该方形数据分析区域的平均标准偏差σ。

4、获取对比度鉴别函数(CDF)曲线，随着方块尺寸的增加，即选取不同的n，建立方块尺寸和平均标准偏差的关系曲线。将不同n下的平均标准偏差表示为占该方形数据分析区域的总的像素的灰度平均值的百分比，建立不同的方块尺寸与该百分比之间的关系曲线即对比度鉴别函数(CDF)曲线。n值变大，即方块尺寸变小，则标准偏差变大，则CDF值变大。如图6所示，横坐标为方块尺寸，竖坐标为相应的CDF值。

Claims (8)

1.一种基于全局二值化的CT图像分辨率自动测试方法，包括采用圆盘标准试件进行测试，其特征在于：包括以下步骤：

1)获取圆盘标准试件的CT图像；

2)对CT图像进行自动阈值分割的全局二值化处理，并且在全局二值化处理后的CT图像中找圆盘标准试件的圆心以及圆盘轮廓线，并且计算圆盘的半径。

3)根据步骤2)中找出的圆盘标准试件的圆心以及圆盘轮廓线和半径，计算CT图像的密度分辨率和空间分辨率。

2.如权利要求1所述的基于全局二值化的CT图像分辨率自动测试方法，其特征在于：上述步骤2)中的自动阈值分割的全局二值化处理包括将步骤1)中测得的图像的灰度值分为0到m级，灰度值为i的像素数为n_i，则总的像素数为各灰度出现的概率为P_i＝n_i/N，将灰度值分为c₀＝{0-k}和c₁＝{k+1-m}，则c₀出现的概率为c₁出现的概率为可以推出c₀平均值为 $u_0=\underset{i=0}{\overset{k}{\Σ}}{\mathrm{iP}}_i/{\mathrm{\ω}}_0=\frac{u\left(k\right)}{\mathrm{\ω}\left(k\right)},$ c₁的平均值为 $u_1=\underset{i=k+1}{\overset{m}{\Σ}}{\mathrm{iP}}_i/{\mathrm{\ω}}_1=\frac{u-u\left(k\right)}{1-\mathrm{\ω}\left(k\right)},$ 其中 $u_k=\underset{i=0}{\overset{k}{\Σ}}{\mathrm{iP}}_i$ 是阈值为k时的灰度平均值，是整体图像的灰度平均值，这样两极之间的方差就可以用下式表示：

s²(k)＝w₀(u-u₀)²+w₁(u₁-u)²

从0到m之间改变k值，计算上式求其最大值便是所求的阈值，然后根据阈值将CT图像二值化，背景为0，圆盘为1；

采用从图像四周向图像中心设置搜索线，查找最外侧的灰度值为1的边界点，作为圆盘的轮廓线；对于获得的圆盘轮廓坐标进行圆周拟合得出圆心，并且计算得出半径。

3.如权利要求2所述的基于全局二值化的CT图像分辨率自动测试方法，其特征在于：上述步骤3)中测试空间分辨率的方法具体包括先获取边界响应函数(ERF)曲线，然后根据先获取边界响应函数(ERF)曲线获得点扩散函数(PSF)曲线，最后根据点扩散函数(PSF)曲线，获得调制传递函数(MTF)曲线。

4.如权利要求3所述的基于全局二值化的CT图像分辨率自动测试方法，其特征在于：获取边界响应函数(ERF)曲线的步骤，具体包括根据步骤2)中确定的圆盘半径R和圆心，选取圆环区域，圆环区域的外径为R+15个像素，内径为R-15个像素，该圆环区域覆盖圆盘标准试件的CT图像的边界，计算圆环区域内所有像素点到圆心的距离，并将距离相同的像素点归为一组，计算每组像素点的灰度值的而平均值，建立每组像素点到圆心的距离与像素灰度值的平均值之间的关系曲线，即边界响应函数(ERF)曲线，并对该边界响应函数(ERF)曲线采用分段插值法和最小二乘立方拟合，获得平滑边界响应函数(ERF)曲线。

5.如权利要求4所述的基于全局二值化的CT图像分辨率自动测试方法，其特征在于：获取点扩散函数(PSF)曲线的具体步骤包括，对获取的边界响应函数(ERF)曲线进行分段拟合，并对每一次拟合得出的多项式进行求导，计算中间点在每个导数解析式对应的值，得出距离与到数值的关系函数，对函数进行归一化处理，得出点扩散函数(PSF)曲线。

6.如权利要求5所述的基于全局二值化的CT图像分辨率自动测试方法，其特征在于：该获取调制传递函数(MTF)曲线为将点扩散函数进行离散傅里叶变换后的幅值在零频处归一化得到。

7.如权利要求2所述的基于全局二值化的CT图像分辨率自动测试方法，其特征在于：步骤3)中测试密度分辨率的方法具体包括，根据步骤2)中获取的圆盘标准试件的圆心和半径，选取图像中心的方形数据分析区域，并且将方形数据分析区域进行n×n分割成方块，n可以为任意自然数，计算每个方块的像素灰度平均值然后通过每个方块内的像素灰度平均值计算整个ROI区域的标准偏差，即得出该n数下该方形数据分析区域的平均标准偏差σ，选取不同的n，将不同n下的平均标准偏差表示为占该方形数据区域的像素的灰度平均值的百分比，建立不同的方块尺寸与该百分比之间的关系曲线即对比度鉴别函数(CDF)曲线。

8.如权利要求1-7中任一项所述的基于全局二值化的CT图像分辨率自动测试方法，其特征在于：该圆盘标准试件的直径选择为工业CT系统成像范围尺寸的70％，厚度不小于20毫米。