一种CT机几何参数校正方法 【技术领域】
本发明属于医学影像处理的技术领域。特别是涉及到一种对x-射线计算机层析成像机几何参数的校正方法。
背景技术
X-射线计算机层析成像机(以下简称CT机)现已成为临床诊断疾病的重要手段，高质量的CT图像可以确认病变细胞的大小和病变细胞的位置。但如果CT机出现图像扭曲，空间分辨率偏低，图像大小不能与实际情况匹配便常常会造成误诊。导致CT图像质量下降的因素很多，几何参数不准就是其中的一个重要因素。
对于第三代CT机而言，旋转半径和中心通道是两个重要的几何参数。从球管中发出x射线束穿透放置在扫描区域中的被检测客体，其衰减值由放置在球管对面的检测器检测。球管发出射线束中穿过旋转中心的一条射线其所对应的检测器其通道号即称为中心通道，而球管到达旋转中心的距离则成为旋转半径。旋转半径不准会导致图像比例失真。而如果中心通道不准会使对面通道无法和实际匹配，这会使图像空间分辨率下降，造成图像扭曲。
目前校正几何参数的方法很少，已知地几种方法可分为直接测量和间接测量两类。直接测量难于实施；间接测量方法通过将需要计算的几何参数转化为测量使用的模体数据，通过实际测量模体数据利用迭代收敛计算几何参数。算法基于实际扫描区域和重建图像之间的对应关系，设计一个目标函数，通过优化算法即当计算的几何参数使目标函数落于极小点时，认为此组参数就是CT机的几何参数。上述算法不仅计算比较费时，而且由于目标函数可能存在不只一个极小点，即使计算得出的几何参数可以满足函数收敛精度，但对图像的重建来说其结果仍不能很好的满足要求。
【发明内容】
本发明提出了一种CT机几何参数的校正方法，其特征在于对针模体进行扫描，利用扫描针模体所得到的生数据坐标中含有的旋转采样和通道角的信息，精确计算出CT机的中心通道和旋转半径，并用所获得的数据校正CT机的几何参数。
本发明所提出的计算中心通道的方法和步骤如下：
(1)扫描单针模体获得正弦图矩阵；
(2)设定阈值，比较正弦图和阈值，获得正弦线；
(3)将通道坐标存入集合DI；
(4)统计DI含有元素的个数SP；
(5)计算中心通道，
本发明所提出的计算旋转半径的方法和步骤如下：
(1)扫描双针模体，获得正弦图矩阵1、2；
(2)测量两个模体中心之间的距离AB；
(3)设定阈值，比较正弦图和阈值，获得两套模型；
(4)调整正弦图矩阵，使扫描的两针模型分别落入第三和第四象限内；
(5)计算两条正弦线中4个特殊点的纵横坐标值；
(6)计算旋转半径，
SO=AB*sin(CI[I1]CNum*CAngel+VI[I3]VNum*2π)sin(VI[I4]-VI[I3]VNum*2π)sin(VI[I4]-VI[I2]VNum*2π+CI[I2]CNum*CAngel)sin(CI[I1]CNum*CAngel)]]>
式中4个特殊点I1、I2、I3和I4为：
I1是正弦线B与纵坐标View＝0的交点；
I2是正弦线A和正弦线B的第一个交点；
I3是正弦线B和中心通道的交点；
I4是正弦线A和中心通道的交点。
在使用时，将本发明的方法对CT机进行测量，并将通过计算所获得的中心通道和旋转半径的数据输入到原CT机的校正表中，则CT机在应用时便会自动按新的更为精确的几何参数对图像进行重建。这样可使图像扭曲程度得到消减，最终使CT图像的质量得到提高。
下面再对本发明的方法作更为详细的论述
1.一、计算中心通道：扫描针模型，将针模型尽量摆放到远离中心的位置，但是不要超过SOV视野。
2.提取针模型在正弦图中的正弦线边界。
所谓正弦线边界提取按照下列步骤实现：
对生数据Raw0进行空气校正得到Raw1，
Raw1=-logRaw0Air]]>
其中，Raw1为灵敏度校正后数据，Air是相同扫描协议下的空气数据
对生数据Raw1进行硬化校正得到Raw2
Raw2=Σi=03hardi(Raw1)i]]>
Raw2为硬化校正后数据，Hardi是机器自带硬化校正表中数据。
对校正结果采用重建的卷积核进行卷积操作得到Raw3，
Raw3＝Raw2*ConvKernel
其中d是采样间隔，n是采样序号
卷积后对每一个View，计算Raw3的差分图像，
计算差分图像Raw3′
Raw3′(i，j)＝Raw3(i，j+1)-Raw3(i，j)其中，i＝0，..，VIEW，j＝0，...，CHANNEL-1，
i，j分别代表旋转角和通道信息。
设定阈值0.5，将差分后的图像Raw3′的阶跃点之间点作为针模型采样点，建立起针模型的采样点集合PIN，完成正弦线边界提取。
3.统计PIN中有元素的数，计为PINNUM.
4.计算中心通道MC
MC=Σi=1PINNUMCI[PIN(i)]PINNUM]]>
其中CI[]表示取括号元素中的通道坐标信息。二、计算旋转半径
为了计算球管到旋转中心的距离即旋转半径SO，采用双针模型。为了精确测量两针之间距离，双针模型可以放置在一个含有两个放针槽，且放针槽可以精确测量的机械装置中。
1、在CT扫描区域中，扫描一个断层平面为圆形的模型，获得正弦图矩阵1；
2、在不同的位置再扫描断层平面为圆形的模型，获得正弦图矩阵2；
3、测量放针槽之间的距离AB，AB距离可通过模体直接精确测量得到；
4、设定阈值，提取正弦线获得两套正弦线矩阵。
5、合并正弦线并调整正弦图矩阵，使扫描的两针模型分别落于第三和第四象限；如图5所示，调整方法如下，寻找纵坐标V0，此V0满足，在V0这行，两条正弦线分别位于中心通道的两侧，同时满足两条正弦线的第一个交点在中心通道右侧；将V0上面的矩阵移到V0下方，定义第一个View在中心通道左侧的正弦线为A，右侧的正弦线为B；
6、计算两条正弦线中4个特殊点I1～I4的横纵坐标值：
I1：正弦线B与纵坐标View＝0的交点，
I2：正弦线A和正弦线B的第一个交点，
I3：正弦线B和中心通道的交点，
I4：正弦线A和中心通道的交点；
7、计算旋转半径
旋转半径SO可以通过下式计算获得。
SO=AB*sin(CI[I1]CNum*CAngel+VI[I3]VNum*2π)sin(VI[I4]-VI[I3]VNum*2π)sin(VI[I4]-VI[I2]VNum*2π+CI[I2]CNum*CAngel)sin(CI[I1]CNum*CAngel)]]>
其中，VI[]表示取括号中点的旋转角度信息也就是纵坐标，CI[]表示取通道角度信息也就是生数据横坐标数值。VNum表示CT旋转一周的采样数，CNum表示CT机一排检测器含有通道的总数，CAngel是扇束张开角
详细的数学推导过程如下：
从图(2)中可知，
SO=ABsin∠AOBsin∠OABsin∠OSBsin(π-∠OSB-∠BOS)]]>
AB是通过摆放得盛放针的机械装置获得，生数据的纵坐标反映了旋转角度信息，横坐标反映了通道角信息。计算角度时，以旋转中心为定点的角度可以利用正弦图的纵坐标获得，以发射源为定点的角度利用正弦图的横坐标获得。例如，计算∠S′OS，当X射线源旋转到S′时，扇形射线束同时穿过A、B两针，也就是说，从正弦图的角度来看，两条正弦线应该相交，因此将交点的纵坐标转换成角度就可以算得角∠S′OS。同样计算角度∠OS′B可以利用交点的横坐标信息。
∠S′OS=VI(I2)VNum*2π]]>
∠OS′B=CI(I2)CNum*CAngel]]>
分别计算公式中的对应的角度，将这些角度转化为可以计算的I1，I2，I3，I4的横纵坐标信息
计算角AOB
∠AOB＝∠AOS-∠BOS
∠AOS=VI(I4)VNum*2π]]>
∠BOS=VI(I3)VNum*2π]]>
∠AOB=VI(I4)-VI(I3)VNum*2π]]>
计算角BAO
∠BAO＝π-∠AOB-∠OBA＝π-∠AOB-(∠BOS-∠BCO)
∠BCO=∠S′OS-∠OS′B=VI(I2)VNum*2π-CI(I2)CNum*CAngel]]>
∠S′OS=VI(I2)VNum*2π]]>
∠OS′B=CI(I2)CNum*CAngel]]>
∠BAO=π-VI(I4)-VI(I3)VNum*2π-(VI(I3)VNum*2π-VI(I2)VNum*2π+CI(I2)CNum*CAngel)]]>
∠&Bgr;&Agr;Ο=πVI(I4)-VI(I2)VNum*2π-CI(I2)CNum*CAngel]]>
计算角OSB
∠OSB=CI(I1)CNum*CAngel]]>
计算角BOS
∠BOS=VI(I3)VNum*2π]]>
利用这种关系，旋转半径可以如下计算
SO=AB*sin(CI[I1]CNum*CAngel+VI[I3]VNum*2π)sin(VI[I4]-VI[I3]VNum*2π)sin(VI[I4]-VI[I2]VNum*2π+CI[I2]CNum*CAngel)sin(CI[I1]CNum*CAngel)]]>
【附图说明】
图1为CT机扫描图像原理图
图2为计算中心通道的流程图；
图3为计算旋转半径的流程图；
图4为单针模型扫描图像，(a)为正弦图，(b)为单针模型图；
图5为双针模型扫描图像，(a)为正弦图，(b)为双针模型图；
图6为应用本发明方法计算的几何参数校正前(a)、后(b)的空间分辨率比较图；
【具体实施方式】
下面将结合具体实施方式对本发明的内容作进一步的说明和补充。虽然实施例是在第三代CT机上应用的，但由于第三代和第四代CT涉及的几何参数可以类似推导，因此方法可以很简单的推广到第四代CT机上去。
实施例1
在单层CT机上应用。机器参数VNum＝1024，CNum＝688，CAngel＝46°*π/180，计算中心通道。
首先使用针模型，工作参数为120KV，150mA，2S，7mm，扫描针模型数据获得1024*688的生数据Raw0，对生数据Raw0进行空气校正得到Raw1，
Raw1=-logRaw0Air]]>
其中，Raw1为灵敏度校正后数据，Air是相同扫描协议下的空气数据
对生数据Raw1进行硬化校正得到Raw2
Raw2=Σi=03hardi(Raw1)i]]>
Raw2为硬化校正后数据，Hardi是机器自带硬化校正表中数据。
对校正结果采用重建的卷积核进行卷积操作得到Raw3，
Raw3＝Raw2*ConvKernel
其中d是采样间隔，n是采样序号
卷积后对每一个View，计算Raw3的差分图像，
计算差分图像Raw3′
Raw3′(i，j)＝Raw3(i，j+1)-Raw3(i，j)其中，i＝0，..，VIEW，j＝0，...，CHANNEL-1，
i，j分别代表旋转角和通道信息。
设定阈值0.5，将差分后的图像Raw3′的阶跃点之间点作为针模型采样点，建立起针模型的采样点集合PIN，将卷积后的生数据进行正弦线提取获得采样点集合PIN，统计集合PIN中有元素的数，计为PINNUM.计算中心通道
MC=Σi=1PINNUMCI[PIN(i)]PINNUM]]>
计算旋转半径，在机械固定装置2中放入针，在不同的放针槽中使用120kv，150mA，2S，7mm扫描金属针。在扫描后的数据中按照上述方法提取正弦线，如图5(b)所示，寻找纵坐标V0，此V0满足在V0这行，两条正弦线分别位于中心通道的两侧，同时满足两条正弦线的第一个交点在中心通道右侧；将V0上面的矩阵移到V0下方，定义第一个View在中心通道左侧的正弦线为A，右侧的正弦线为B；计算两条正弦线中4个特殊点I1～I4的横纵坐标值：I1是正弦线B与纵坐标View＝0的交点；I2是正弦线A和正弦线B的第一个交点；I3是正弦线B和中心通道的交点；I4是正弦线A和中心通道的交点。进行矩阵的平移操作(1024*688)计算CI[I1]，VI[I3]，VI[I4]VI[I2]CI[I2]
计算SO=AB*sin(CI[I1]CNum*CAngel+VI[I3]VNum*2π)sin(VI[I4]-VI[I3]VNum*2π)sin(VI[I4]-VI[I2]VNum*2π+CI[I2]CNum*CAngel)sin(CI[I1]CNum*CAngel)]]>
RawI(i,j)=-logRaw0(i,j)I0]]>其中I0是仅扫描空气的数据
RawI(i,j)=Σk=04ak(Raw1)k,]]>其中ai是硬化校正表
利用上述方法进行计算中心通道和旋转半径。结果显示CT扇束张角为46°，扫描旋转一周采样1024次，每次采样有688个检测器通道进行检测，扫描电压120kv，切片厚度7mm，扫描电流150mA，扫描时间2s。中心通道原值343.32，校正后值343.27。旋转半径原值593.7，校正后值为592.6，将校正后的这组几何参数输入CT机的校正表中，便可使用。
实施例2
在双层CT机上应用。机器参数VNum＝1024，CNum＝720，CAngel＝46°*π/180，计算中心通道。
首先使用针模型，工作参数为140KV，100mA，1S，5mm，扫描针模型数据获得1024*720的生数据Raw0，对生数据Raw0进行空气校正得到Raw1，
Raw1=-logRaw0Air]]>
其中，Raw1为灵敏度校正后数据，Air是相同扫描协议下的空气数据
对生数据Raw1进行硬化校正得到Raw2
Raw2=Σi=03hardi(Raw1)i]]>
Raw2为硬化校正后数据，Hardi是机器自带硬化校正表中数据。
对校正结果采用重建的卷积核进行卷积操作得到Raw3，
Raw3＝Raw2*ConvKernel
其中d是采样间隔，n是采样序号卷积后对每一个View，计算Raw3的差分图像，
计算差分图像Raw3′
Raw3′(i，j)＝Raw3(i，j+1)-Raw3(i，j)其中，i＝0，..，VIEW，j＝0，...，CHANNEL-1，
i，j分别代表旋转角和通道信息。
设定阈值0.2，将差分后的图像Raw3′的阶跃点之间点作为针模型采样点，建立起针模型的采样点集合PIN，将卷积后的生数据进行正弦线提取获得采样点集合PIN，统计集合PIN中有元素的数，计为PINNUM.计算中心通道
MC=Σi=1PINNUMCI[PIN(i)]PINNUM]]>
计算旋转半径，在机械固定装置2中放入针，在不同的放针槽中使用140kv，100mA，1S，5mm扫描金属针。在扫描后的数据中按照上述方法提取正弦线，如图5(b) 所示，寻找纵坐标V0，此V0满足在V0这行，两条正弦线分别位于中心通道的两侧，同时满足两条正弦线的第一个交点在中心通道右侧；将V0上面的矩阵移到V0下方，定义第一个View在中心通道左侧的正弦线为A，右侧的正弦线为B；计算两条正弦线中4个特殊点I1～I4的横纵坐标值：I1：正弦线B与纵坐标View＝0的交点；I2：正弦线A和正弦线B的第一个交点；I3：正弦线B和中心通道的交点；I4：正弦线A和中心通道的交点。进行矩阵的平移操作(1024*688)计算CI[I1]，VI[I3]，VI[I4]VI[I2]CI[I2]
计算SO=AB*sin(CI[I1]CNum*CAngel+VI[I3]VNum*2π)sin(VI[I4]-VI[I3]VNum*2π)sin(VI[I4]-VI[I2]VNum*2π+CI[I2]CNum*CAngel)sin(CI[I1]CNum*CAngel)]]>
RawI(i,j)=-logRaw0(i,j)I0]]>其中I0是仅扫描空气的数据
RawI(i,j)=Σk=04ak(Raw1)k,]]>其中ai是硬化校正表
利用上述方法进行计算中心通道和旋转半径。结果显示CT扇束张角为46°，扫描旋转一周采样1024次，每次采样有720个检测器通道进行检测，扫描电压140kv，切片厚度1mm，扫描电流100mA，扫描时间1s。中心通道原值359.41，校正后值359.24。旋转半径原值616.3，校正后值为614.7，将校正后的这组几何参数输入CT机的校正表中，便可使用。
采用本发明的方法校正后的几何参数，在同一CT机上，使用相同扫描条件120KV，150mA，2s，7mm参数扫描AAPM空间分辨率模体。空间分辨率模体如图6所示，其中图6a为应用原来系统采用的几何参数重建结果，图6b为应用本方法计算得到的几何参数重建结果，结果显示前者空间分辨率能分辨到第五行0.6mm，而后者空间分辨率可以显示第六行，0.5mm。同时图6中a、b两幅图像比较也可以看出，图6b点圆度优于图6a，新参数消减了图像的扭曲。从结果来看，图像的空间分辨率得到了提高，图像扭曲得到消减，经过几何参数校正后的图像质量得到了提高。