一种改进圆盘卡法的CT空间分辨率测量方法与流程

本发明属于图像处理
技术领域：
，涉及一种改进圆盘卡法的ct(computedtomography)空间分辨率测量方法。
背景技术：
：空间分辨率又称为高对比度分辨率，它是描述ct系统的重要指标，目前的测试方法有线对卡法、圆孔卡法和圆盘卡法。线对卡法和圆孔卡法虽然测试过程简单，但对制作要求高。圆盘卡制作简单，但测试过程复杂，测试结果需要繁杂的计算。圆盘卡法的计算过程中由于噪声和圆盘卡制作精度不足，这对于在测试过程中最为重要的一步即对于平均灰度计算的影响较大，最终得到的空间分辨率的结果会与真实值有一定的差距。为了减少这种影响，本发明改进了圆盘卡法的测试过程，使其空间分辨率更加准确，更能反映ct系统的性能，因此对提高圆盘卡在空间分辨率测量方法精度的研究具有较大的实际意义。技术实现要素：有鉴于此，本发明的目的在于提供一种改进圆盘卡法的ct空间分辨率测量方法，用于在计算圆盘卡边缘灰度值时降低由噪声和圆盘卡制作精度不足对边缘平均灰度值计算的影响，从而提高在计算工业ct系统的空间分辨率时的准确度。为达到上述目的，本发明提供如下技术方案：一种改进圆盘卡法的ct空间分辨率测量方法，具体包括以下步骤：s1：安装并启动ct检测装置，使得射线源(1)产生的有合适的角度的扇形束射线能够覆盖圆盘卡的全部区域，这个角度设为2γ；s2：扇形束射线扫描，使射线源(1)与曲面探测器(2)绕圆盘卡(3)的圆心旋转180°+2γ来获得完整的投影数据，并传送到控制及图像处理系统(4)中存储；s3：重建圆盘卡的二维图像：根据得到的投影数据利用扇形束fbp算法重建圆盘卡的二维图像；s4：利用rtv算法来得到去噪后的灰度图像；s5：使用改进的圆盘卡法对重建的二维图像和去噪后的灰度图像进行处理；s6：得到圆盘卡的mtf曲线和ct空间分辨率。进一步：步骤s1中所述的检测装置包括射线源(1)、曲面探测器(2)以及控制及图像处理系统(4)，所述射线源(1)、曲面探测器(2)的信号线路与控制与图像处理系统(4)相连，射线源(1)中心与圆盘卡(3)的圆心连线与曲面探测器(2)大致垂直，射线源(1)与曲面探测器(2)绕圆盘卡(3)圆心旋转，使得射线源(1)产生的有合适的角度的扇形束射线能够覆盖圆盘卡的全部区域，这个角度设为2γ。进一步，所述步骤s2具体包括：在控制及图像处理系统(4)的控制下，首先将射线源(1)中心对准曲面探测器(2)的中心，射线源(1)与圆盘卡(3)的圆心连线与曲面探测器(2)大致垂直，射线源与曲面探测器(2)绕圆盘卡(3)的圆心至少旋转180°+2γ来获得完整的投影数据，然后传送到控制及图像处理系统(4)中存储。进一步，步骤s3中，所述利用等距扇形束fbp算法重建圆盘卡的二维图像的计算公式为：其中，(x,y)表示待重建点的二维坐标，f(x,y)为经过计算得到的坐标为(x,y)的重建值，r(γ,β)为射线源与旋转中心连线与纵轴之间的夹角为β且扇角为γ时的投影值，d为射线源到旋转中心的距离，d′为重建点与射线源之间的距离，γ′为射线源与旋转中心连线和重建点与射线源连线的夹角，γm为扇形束的最大幅角，即扇形束角度的一半，w为频率变量。进一步，步骤s4中，所述利用rtv算法来得到去噪后的灰度图像的计算公式为：其中，i表示输入图像，p表示2d图像像素点的索引，s表示输出结构图像，ε固定为0.001；λ表示一个不可或缺的权重，用来控制图像的光滑程度，但是仅仅调节它不会使纹理分离太多，而增加λ也会造成图像的模糊并且纹理反而保留下来，一般λ选取为0.01到0.03之间；及其中，q表示以p点为中心的一个正方形区域r(p)内所有的像素点的索引，r(p)为以p点为中心给定长度的正方形区域，此处长度根据实际情况设置；和分别表示输出结构图像s在待重建点的二维坐标x和y方向的偏导在q点处的值；高斯核函数：其中，xp、xq和yp、yq分别为p点和q点在图像中的横纵坐标，σ为r(p)区域灰度图像的标准方差，exp为自然数e。进一步，所述步骤s5具体包括以下步骤：s51：采用i表示输入图像，i(i,j)代表图像i中坐标为(i,j)的位置的灰度；s52：将输入图像i进行归一化处理：i＝(i-min(i(:)))/(max(i(:))-min(i(:)))，绘制出灰度图像中线的灰度变化图，选取合适阈值，进行阈值分割，分割的圆盘卡部分灰度设定为1，其它部位为0；s53：计算圆盘卡的半径，采用形态学腐蚀的方法：首先令ir＝i，k＝0(该步不参与循环)，然后循环下述步骤：ir1(i,j)＝min{[ir(i-1,j),ir(i+1,j),ir(i,j),ir(i,j+1),ir(i,j-1)])}；ir＝ir1；k＝k+1；每循环一次判断直至i＝0，圆盘卡的半径r＝k-1；ir1(i,j)表示图像ir1中横坐标为i，纵坐标为j的点的灰度值；s54：计算圆盘卡边缘轮廓的平均灰度：选取合适的t，即计算圆盘卡边缘向内和向外t层，每层的平均灰度值；计算过程如下：首先令io＝i(该步不参与循环)；i1(i,j)＝min([i(i-1,j),i(i+1,j),i(i,j),i(i,j+1),i(i,j-1)])；io1(i,j)＝max([io(i-1,j),io(i+1,j),io(i,j),io(i,j+1),io(i,j-1)])；判断i1与i中灰度值不同的像素点的位置，计算这些位置对应在原图(经过fbp重建得到的，未经过rtv去噪的图像)像素点的平均灰度值作为向内腐蚀的第一层环的平均灰度，同理判断io1与io中灰度值不同的像素点的位置，计算这些位置对应在原图像素点的平均灰度值作为向外腐蚀的第1层环的平均灰度；再让i＝i1，io1＝io，重复上面的步骤计算向内和向外第2层的平均灰度，如此重复计算出圆盘卡边缘向内和向外t层，每层的平均灰度；其中i1(i,j)和io1(i,j)分别表示图像i1和io1中横坐标为i，纵坐标为j的点的灰度值。进一步，所述步骤s6具体包括以下步骤：s61：绘制圆盘边缘轮廓的灰度变化图，横坐标为半径从小到大，与s53计算的圆盘半径和s54中t的取值有关，范围为r-t+1到r+t；纵坐标为对应的平均灰度，由s54计算得到；s62：得到边缘响应函数erf；具体包括：按照推荐的拟合点数，依次选取相应的点数组合；第二个组合的第一个点是第一个组合的第二个点；对每个组合做最小二乘拟合，用拟合所得的中间点代替原来点的灰度值，依次重复操作，计算出全部拟合后的灰度值，从而得到距离和拟合灰度值之间的关系；删除开始端和结束端多余数据，根据距离和拟合灰度值之间的关系得到erf；s63：对erf生成的结果再做分段拟合，类似s62，并对每一组拟合得到的多项式求导，再由每个导数解析式计算中间点的导数值，得到距离和导数值之间的关系；s64：用最大导数值对所有导数值进行归一化处理得到点扩展函数psf；s65：得到mtf曲线：计算psf的傅里叶变化，得到傅里叶变化的振幅；对振幅随频域变化曲线在零频率处归一化处理，得到mtf曲线；s66：ct系统的空间分辨率是调制度为10％时对应的线对数，根据s66得到的mtf曲线中直接得到。本发明的有益效果在于：本发明所述方法在计算圆盘卡边缘灰度值时降低了由噪声和圆盘卡制作精度不足的情况下的对边缘平均灰度值计算的影响，从而在计算ct系统的空间分辨率时具有更加准确合理的结果。本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。附图说明为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作优选的详细描述，其中：图1为本发明的工业部件扫描结构示意图；图2为等间角扇形束下重建点c的几何关系图图3为本发明所述ct空间分辨率测量方法的流程图；图4为真实数据圆盘卡重建图像；图5为采用传统圆盘卡法得到的erf曲线图；图6为采用本发明所述方法得到的erf曲线图；附图标记：1-射线源，2-曲面探测器，3-圆盘卡，4-控制及图像处理系统。具体实施方式以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。其中，附图仅用于示例性说明，表示的仅是示意图，而非实物图，不能理解为对本发明的限制。请参阅图1～图3，图1为本发明采用的工业部件扫描结构示意图，如图1所示，该检测装置包括射线源1、曲面探测器2和控制及图像处理系统4。射线源1、曲面探测器2的信号线路与控制及图像处理系统4相连，射线源1中心与圆盘卡3的圆心连线与曲面探测器2大致垂直，射线源1与曲面探测器2绕圆盘卡3圆心旋转，使得射线源1产生的扇形束射线能够覆盖圆盘卡的全部区域。图2为等角度扇形束下重建点c的几何关系图，如图2所示，o为旋转中心s为射线源，β为射线源与旋转中心o连线与纵轴之间的夹角，r为重建点c到旋转中心的距离，为旋转中心与重建点与x轴的夹角。图3为本发明所述基于圆盘卡的空间分辨率测量方法的流程图，如图3所示，该方法具体包括以下步骤：s1：安装检测装置：该检测装置包括射线源1、曲面探测器2以及控制及图像处理系统4，所述射线源1、曲面探测器2的信号线路与控制与图像处理系统4相连，射线源1中心与圆盘卡3的圆心连线与曲面探测器2大致垂直，射线源1与曲面探测器2绕圆盘卡3圆心旋转，使得射线源1产生的有合适的角度的扇形束射线能够覆盖圆盘卡的全部区域，这个角度设为2γ。s2：扇形束射线扫描，获得完整的投影数据：在控制及图像处理系统4的控制下，首先将射线源1中心对准曲面探测器2的中心，射线源1与圆盘卡3的圆心连线与曲面探测器2大致垂直，射线源与曲面探测器2绕圆盘卡3的圆心旋转180°+2γ来获得完整的投影数据，然后传送到控制及图像处理系统4中存储。s3：重建圆盘卡的二维图像：根据得到的投影数据利用扇形束fbp算法重建圆盘卡的二维图像，计算公式为：其中，(x,y)表示待重建点的二维坐标，f(x,y)为经过计算得到的坐标为(x,y)的重建值，r(γ,β)为射线源与旋转中心连线与纵轴之间的夹角为β且扇角为γ时的投影值，d为射线源到旋转中心的距离，d′为重建点与射线源之间的距离，γ′为射线源与旋转中心连线和重建点与射线源连线的夹角，γm为扇形束的最大幅角，即扇形束角度的一半，w为频率变量。s4：利用rtv算法来得到去噪后的灰度图像，计算公式为：其中，i表示输入图像，p表示2d图像像素点的索引，s表示输出结构图像，ε固定为0.001；λ表示一个不可或缺的权重，用来控制图像的光滑程度，但是仅仅调节它不会使纹理分离太多，而增加λ也会造成图像的模糊并且纹理反而保留下来，一般λ选取为0.01到0.03之间；及其中，q表示以p点为中心的一个正方形区域r(p)内所有的像素点的索引，r(p)为以p点为中心给定长度的正方形区域，此处长度根据实际情况设置；和分别表示输出结构图像s在待重建点的二维坐标x和y方向的偏导在q点处的值；高斯核函数：其中，xp、xq和yp、yq分别为p点和q点在图像中的横纵坐标，σ为r(p)区域灰度图像的标准方差，exp为自然数e。s5：使用改进的圆盘卡法对重建的二维图像和去噪后的灰度图像进行处理，具体步骤为：s51：采用i表示输入图像，i(i,j)代表图像i中坐标为(i,j)的位置的灰度；s52：将输入图像i进行归一化处理：i＝(i-min(i(:)))/(max(i(:))-min(i(:)))，绘制出灰度图像中线的灰度变化图，选取合适阈值，进行阈值分割，分割的圆盘卡部分灰度设定为1，其它部位为0；s53：计算圆盘卡的半径，采用形态学腐蚀的方法：首先令ir＝i，k＝0(该步不参与循环)，然后循环下述步骤：ir1(i,j)＝min{[ir(i-1,j),ir(i+1,j),ir(i,j),ir(i,j+1),ir(i,j-1)])}；ir＝ir1；k＝k+1；每循环一次判断直至i＝0，圆盘卡的半径r＝k-1；ir1(i,j)表示图像ir1中横坐标为i，纵坐标为j的点的灰度值；s54：计算圆盘卡边缘轮廓的平均灰度：选取合适的t，一般取100，即计算圆盘卡边缘向内和向外100层，每层的平均灰度值；计算过程如下：首先令io＝i(该步不参与循环)；i1(i,j)＝min([i(i-1,j),i(i+1,j),i(i,j),i(i,j+1),i(i,j-1)])；io1(i,j)＝max([io(i-1,j),io(i+1,j),io(i,j),io(i,j+1),io(i,j-1)])；判断i1与i中灰度值不同的像素点的位置，计算这些位置对应在原图(经过fbp重建得到的，未经过rtv去噪的图像)像素点的平均灰度值作为向内腐蚀的第一层环的平均灰度，同理判断io1与io中灰度值不同的像素点的位置，计算这些位置对应在原图像素点的平均灰度值作为向外腐蚀的第1层环的平均灰度；再让i＝i1，io1＝io，重复上面的步骤计算向内和向外第2层的平均灰度，如此重复计算出圆盘卡边缘向内和向外100层，每层的平均灰度；其中i1(i,j)和io1(i,j)分别表示图像i1和io1中横坐标为i，纵坐标为j的点的灰度值。s6：计算erf和psf曲线，得到圆盘卡的mtf曲线和ct空间分辨率，具体步骤为：s61：绘制圆盘边缘轮廓的灰度变化图，横坐标为半径从小到大，与s53计算的圆盘半径和s54中t的取值有关，范围为r-t+1到r+t；纵坐标为对应的平均灰度，由s54计算得到；s62：得到边缘响应函数erf：按照表1(见说明书附图)中推荐的拟合点数，依次选取相应的点数组合；第二个组合的第一个点是第一个组合的第二个点；对每个组合做最小二乘拟合，用拟合所得的中间点代替原来点的灰度值，依次重复操作，计算出全部拟合后的灰度值，从而得到距离和拟合灰度值之间的关系；删除开始端和结束端多余数据，根据距离和拟合灰度值之间的关系得到erf。表1推荐适用的各项参数(单位：像素)图像尺寸圆盘图像直径方块最大尺度像素距离单位拟合点数256235120.10011512470240.050211024940480.0254120481880960.012581表1为本发明采用的边缘响应函数过程中推荐的拟合点数。根据表1中推荐的拟合点数计算边缘响应函数，不包括表中尺寸的可选相近尺寸的拟合点数。s63：对erf生成的结果再做分段拟合，类似s62，并对每一组拟合得到的多项式求导，再由每个导数解析式计算中间点的导数值，得到距离和导数值之间的关系；s64：用最大导数值对所有导数值进行归一化处理得到点扩展函数psf；s65：得到mtf曲线：计算psf的傅里叶变化，得到傅里叶变化的振幅；图像的截至频率定义为0.5像素/线对，傅里叶变换后的最高频率应不低于图像矩阵截至频率的4倍；按照采样定律，psf的采样间隔不大于0.25像素，为了获得光滑的mtf曲线，频域内的采样间隔应小于0.01lp/pixel(亦即对于psf的采样应大于100个像素)；对振幅随频域变化曲线在零频率处归一化处理，得到mtf曲线；s66：ct(computedtomography，计算机断层成像)系统的空间分辨率是调制度为10％时对应的线对数，根据s59得到的mtf曲线中直接得到。实验对比分析：针对图4所示的真实数据圆盘卡重建图像进行分析，在圆盘卡重建后图像受到较大噪声的影响时，传统的圆盘卡法计算出的erf图如图5，曲线震荡，很不稳定，无法进行后续空间分辨率的计算。而本发明提出的改进的圆盘卡法可以得到较光滑的曲线，如图6所示，曲线变化符合物理逻辑，有利于降低噪声干扰计算出更准确的ct空间分辨率。最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。当前第1页12