一种基于交叉视觉皮质模型的医学图像配准方法，属于基于灰度的图像配准方法，解决现有医学图像配准方法运行时间较长的问题，对现有的交叉视觉皮质模型进行改进，减少其卷积运算，使之具有更快的运行速度。本发明包括边缘检测步骤、提取特征向量步骤、计算配准参数步骤和配准步骤。本发明利用点火次数矩阵的重心特征点实现图像配准，减少了运行时间，同时利用边缘检测将待配准图像和参考图像中的像素分类为边缘像素和非边缘像素，并对边缘像素采用原始ICM公式计算像素状态值，对非边缘像素采用改进的ICM公式计算像素状态值，减少了公式中的卷积运算，从而进一步降低了运行时间。
1.一种基于交叉视觉皮质模型的医学图像配准方法，适用于保形变换，包括下述步骤： 一.边缘检测步骤：
输入待配准图像H和参考图像R，分别利用sobel横向边缘提取算子Gx和sobel纵向边缘提取算子Gy对H和R中每个像素的8邻域进行卷积运算，提取目标边缘，得到待配准边缘图像H1和参考边缘图像R1； 二.提取特征向量步骤，包括下述子步骤：
(2.1)初始化，设置运行次数N＝20～50，置运行次数变量n＝1，各像素状态值Fij[0]＝0，响应矩阵Y[0]各元素值Yij[0]＝0，各像素阈值θij[0]＝0，点火次数矩阵M[0]各元素值Mij[0]＝0，i、j分别为待配准图像H和参考图像R中像素的横坐标和纵坐标； (2.2)计算像素状态值Fij[n]： 对H和R中每个像素计算像素状态值Fij[n]：
判断H和R中的像素是否分别为H1和R1中目标边缘点， 是则：Fij[n]＝f×Fij[n-1]+Sij+Wkl{Ykl[n-1]}， 否则：Fij[n]＝f×Fij[n-1]+Sij+C×Yij[n-1]， 式中，权值矩阵元素值
Sij为该像素灰度值，k、1分别为当前像素8邻域内像素的横坐标和纵坐标，权值常数C为Wk1所有元素之和，0＜状态衰减系数f＜1；Yk1[n]为当前像素的八邻域像素响应矩阵元素值，Yk1[n]∈Y[n]； θij[n]＝g×θij[n-1]+h×Yij[n-1]， 其中，θij[n]为各像素阈值；0＜阈值衰减系数g＜1，10＜阈值常数h＜100； (2.3)计算点火次数矩阵M[n]各元素值Mij[n]： 当Yij[n]＝1时，Mij[n]＝Mij[n-1]+1， 当Yij[n]＝0时，Mij[n]＝Mij[n-1]； 根据上式，定义待配准图像H的点火次数矩阵M[n]为MH[n]，参考图像R[n]的点火次数矩阵M[n]为MR[n]； (2.4)计算重心特征点；
分别对H和R计算MH[n]和MR[n]矩阵的重心特征点PH[n]和PR[n]： PH [n]的x坐标和y坐标分别为：
(2.5)置n＝n+1，判断是否n＞N，是则转子步骤(2.6)，否则转子步骤(2.2)； (2.6)将每次迭代得到的PH[n]和PR[n]分别按n的顺序从上至下依次排列得到N×2的待配准图像重心特征矩阵 和参考图像重心特征矩阵 转步骤三； 三.计算配准参数步骤：
(3.1)计算H和R之间的x坐标平移参数Δx和y坐标平移参数Δy： (3.2)计算旋转参数Δθ：
其中，反正切分子向量 反正切分母向量 向量内元素序号p＝1～N； 四.配准步骤：
式中，x和y为H的像素坐标，x′和y′表示平移变换后待配准图像H′中像素的坐标； (4.2)对H′进行旋转变换：
得到最终配准图像H″，其中x″和y″表示最终配准图像H″中像素的坐标。
2.如权利要求1所述的医学图像配准方法，其特征在于，所述提取特征向量步骤的子步骤(2.2)中： 先选定所述状态衰减系数f，再选定所述阈值衰减系数g，最后选定所述 阈值常数h； 选定所述状态衰减系数f的方式为：首先设定阈值衰减系数g为0.9，设定阈值常数h＝20，让f以步长0.1从0.1变化至0.9，根据g、h和各f值计算对应的9个最终配准图像H″，分别计算每个H″与R之间的均方差MSE： 式中：I，J分别表示H″与R的长和宽；
选取其中最小MSE值对应的最终配准图像H″，其对应的f值为选定的状态衰减系数； 选定所述阈值衰减系数g的方式为：选定状态衰减系数f后，设定阈值常数h＝20，让g以步长0.1从0.9变化至0.1，根据f、h和各g值计算对应的9个最终配准图像H″，分别计算每个H″与R之间的均方差MSE；选取其中最小MSE值对应的最终配准图像H″，其对应的g值为选定的阈值衰减系数； 选定所述阈值常数h的方式为：选定状态衰减系数f和阈值衰减系数g后，让h以步长
10从10变化至100，根据f、g和各h值计算对应的10个最终配准图像H″，分别计算每个H″与R之间的均方差MSE；选取其中最小MSE值对应的最终配准图像H″，其对应的h值为选定的阈值常数。
一种基于交叉视觉皮质模型的医学图像配准方法
技术领域
[0001] 本发明属于基于灰度的图像配准方法，具体涉及一种基于交叉视觉皮质模型的医学图像配准方法。
背景技术
[0002] 图像配准技术在遥感数据分析、计算机视觉、医学图像处理等领域都有很广泛的应用，是目标识别、图像融合、时序图像分析、变化检测等实际应用的关键步骤，在环境检测、天气预报等应用领域都有不可替代的地位。
[0003] 根据所采用的图像信息，现有图像配准方法可以分为两大类：基于特征的图像配准方法和基于灰度的图像配准方法。
[0004] 基于特征的图像配准方法，例如：中国专利号200680048083.5，名称为“基于点的自适应弹性图像配准”、中国专利号200810031575.4，名称为“基于直线特征图像配准中的特征匹配方法”。这类图像配准方法，依赖于特征的提取，以图像中不变特征，如点(包括角点、高曲率点等)、线、边缘、轮廓、闭合区域以及统计特征不变量比如重心等，作为图像配准的特征，要求特征提取具有可靠性和鲁棒性。
[0005] 基于灰度的图像配准方法，例如：中国专利号200810019451.4，名称为“基于量子行为粒子群算法的多分辨率医学图像配准方法”，公布了利用归一化互信息作为目标函数的配准方法；中国专利号200710052491.4，名称为“一种多相似性测度图像配准方法”。这类图像配准方法，与图像的像素灰度值密切相关，不用对图像进行特征提取，通常以一定的目标函数作为测度，通过优化方法寻找最优情况下的配准参数，运行时间会有所增加。
[0006] 交叉视觉皮质模型(Intersecting Cortical Model，ICM)是在脉冲耦合神经网络模型(Pulse Coupled Neural Network，PCNN)上的简化和改进，这两者被称之为第三代神经网络模型，已成功地被应用于图像分割、图像去噪、特征提取等领域。ICM具有的动态阈值、同步脉冲发放等特征非常适合于图像处理，而且同传统神经网络模型相比，它具有自适应、不需要训练和学习的优点。同时，ICM还具有平移、旋转不变性等特点。ICM模型的数学方程为：
[0007] Fij[n]＝f×Fij[n-1]+Sij+Wkl{Ykl[n-1]}，
[0008]
[0009] θij[n]＝g×θij[n-1]+h×Yij[n-1]，
[0010] 式中，权值矩阵元素值
[0011] Fij[n]为像素状态值，Sij为像素灰度值，i、j分别为图像中像素的横坐标和纵坐标，k、1分别为当前像素8邻域内像素的横坐标和纵坐标，f为状态衰减系数，Yk1[n]为当前像素的8邻域像素响应矩阵元素值，Yk1[n]∈Y[n]；θij[n]为各像素阈值，g为阈值衰减系数，h为阈值常数。
[0012] 中国专利号200910086060.9，名称为“一种基于改进交叉视觉皮质模型的图像分割方法”，公布了交叉视觉皮质模型在图像分割中的应用。虽然ICM在图像处理已经有广泛的应用，但在图像配准领域还没有出现过相关应用。图像配准过程中有大量的重复的乘法和加法运算，ICM本身也存在卷积运算，导致运行时间较长。
[0013] 本发明中，保形变换是指待配准图像与参考图像之间只存在平移，旋转以及图像大小变化的变换。
发明内容
[0014] 本发明提供一种基于交叉视觉皮质模型的医学图像配准方法，解决现有医学图像配准方法运行时间较长的问题，对现有的交叉视觉皮质模型进行改进，减少其卷积运算，使之具有更快的运行速度。
[0015] 本发明所提供的一种基于交叉视觉皮质模型的医学图像配准方法，适用于保形变换，包括下述步骤：
[0016] 二.边缘检测步骤：
[0017] 输入待配准图像H和参考图像R，分别利用sobel横向边缘提取算子Gx和sobel纵向边缘提取算子Gy对H和R中每个像素的8邻域进行卷积运算，提取目标边缘，得到待配准边缘图像H1和参考边缘图像R1；
[0018]
[0019] 二.提取特征向量步骤，包括下述子步骤：
[0020] (2.1)初始化，设置运行次数N＝20～50，置运行次数变量n＝1，各像素状态值Fij[0]＝0，响应矩阵Y[0]各元素值Yij[0]＝0，各像素阈值θij[0]＝0，点火次数矩阵M[0]各元素值Mij[0]＝0，i、j分别为待配准图像H和参考图像R中像素的横坐标和纵坐标；
[0021] (2.2)计算像素状态值Fij[n]：
[0022] 对H和R中每个像素计算像素状态值Fij[n]：
[0023] 判断H和R中的像素是否分别为H1和R1中目标边缘点，
[0024] 是则：Fij[n]＝f×Fij[n-1]+Sij+Wkl{Ykl[n-1]}，
[0025] 否则：Fij[n]＝f×Fij[n-1]+Sij+C×Yij[n-1]，
[0026] 式中，权值矩阵元素值
[0027] Sij为该像素灰度值，k、1分别为当前像素8邻域内像素的横坐标和纵坐标，权值常数C为Wk1所有元素之和，0＜状态衰减系数f＜1；Yk1[n]为当前像素的八邻域像素响应矩阵元素值，Yk1[n]∈Y[n]；
[0028]
[0029] θij[n]＝g×θij[n-1]+h×Yij[n-1]，
[0030] 其中，θij[n]为各像素阈值；0＜阈值衰减系数g＜1，10＜阈值常数h＜100；
[0031] (2.3)计算点火次数矩阵M[n]各元素值Mij[n]：
[0032] 当Yij[n]＝1时，Mij[n]＝Mij[n-1]+1，
[0033] 当Yij[n]＝0时，Mij[n]＝Mij[n-1]；
[0034] 根据上式，定义待配准图像H的点火次数矩阵M[n]为MH[n]，参考图像R[n]的点火次数矩阵M[n]为MR[n]；
[0035] (2.4)计算重心特征点；
[0036] 分别对H和R计算MH[n]和MR[n]矩阵的重心特征点PH[n]和PR[n]：
[0037] PH[n]的x坐标和y坐标分别为：
[0038]
[0039] PR[n]的x坐标和y坐标分别为：
[0040]
[0041] (2.5)置n＝n+1，判断是否n＞N，是则转子步骤(2.6)，否则转子步骤(2.2)；
[0042] (2.6)将每次迭代得到的PH[n]和PR[n]分别按n的顺序从上至下依次排列得到N×2的待配准图像重心特征矩阵 和参考图像重心特征矩阵 转步骤三；
[0043] 三.计算配准参数步骤：
[0044] (3.1)计算H和R之间的x坐标平移参数Δx和y坐标平移参数Δy：
[0045]
[0046] (3.2)计算旋转参数Δθ：
[0047]
[0048]
[0049]
[0050] 其中，反正切分子向量 反正切分母向量 向量内元素序号p＝1～N；
[0051] 四.配准步骤：
[0052] (4.1)对待配准图像H进行平移变换：
[0053]
[0054] 式中，x和y为H的像素坐标，x′和y′表示平移变换后待配准图像H′中像素的坐标；
[0055] (4.2)对H′进行旋转变换：
[0056]
[0057] 得到最终配准图像H″，其中x″和y″表示最终配准图像H″中像素的坐标。
[0058] 所述的医学图像配准方法，其特征在于，所述提取特征向量步骤的子步骤(2.2)中：
[0059] 先选定所述状态衰减系数f，再选定所述阈值衰减系数g，最后选定所述阈值常数h；
[0060] 选定所述状态衰减系数f的方式为：首先设定阈值衰减系数g为0.9，设定阈值常数h＝20，让f以步长0.1从0.1变化至0.9，根据g、h和各f值计算对应的9个最终配准图像H″，分别计算每个H″与R之间的均方差MSE：
[0061]
[0062] 式中：I，J分别表示H″与R的长和宽；
[0063] 选取其中最小MSE值对应的最终配准图像H″，其对应的f值为选定的状态衰减系数；
[0064] 选定所述阈值衰减系数g的方式为：选定状态衰减系数f后，设定阈值常数h＝20，让g以步长0.1从0.9变化至0.1，根据f、h和各g值计算对应的9个最终配准图像H″，分别计算每个H″与R之间的均方差MSE；选取其中最小MSE值对应的最终配准图像H″，其对应的g值为选定的阈值衰减系数；
[0065] 选定所述阈值常数h的方式为：选定状态衰减系数f和阈值衰减系数g后，让h以步长10从10变化至100，根据f、g和各h值计算对应的10个最终配准图像H″，分别计算每个H″与R之间的均方差MSE；选取其中最小MSE值对应的最终配准图像H″，其对应的h值为选定的阈值常数。
[0066] 本发明的主要优点在于利用边缘检测将待配准图像和参考图像中的像素分类为边缘像素和非边缘像素，并对边缘像素采用原始ICM公式计算像素状态值，对非边缘像素采用改进的ICM公式计算像素状态值，其中改进的ICM公式减少了公式中的卷积运算，从而降低了公式的运行时间，同时还提出了利用点火次数矩阵的重心特征点实现图像配准的思路，有效的利用了神经元本身特性。表1为采用粒子群(PSO)方法、鲍威尔(Powell)方法、仅采用点火次数矩阵M的配准方法以及本发明在配准结果上的对比。从表中可以看到仅采用点火次数矩阵M的配准方法以及本发明在配准结果中准确度高于粒子群(PSO)方法和鲍威尔(Powell)方法的配准结果。
[0067] 表1
[0068]
[0069] 表1中，第一列变换参数括号内的参数依次为：参考图像R和待配准图像H之间理论横坐标平移量、理论纵坐标平移量、理论旋转角度；第二列粒子群(PSO)方法配准结果括号内的参数依次为：待配准图像H和采用PSO配准方法得到的最终配准图像H″之间实际横坐标平移量、实际纵坐标平移量、实际旋转角度；第三列鲍威尔(Powell)方法配准结果括号内的参数依次为：待配准图像H和采用鲍威尔(Powell)配准方法得到的最终配准图像H″之间实际横坐标平移量、实际纵坐标平移量、实际旋转角度；第四列仅采用点火次数矩阵M的配准方法配准结果括号内的参数依次为：待配准图像H和采用点火矩阵配准方法得到的最终配准图像H″之间实际横坐标平移量、实际纵坐标平移量、实际旋转角度；第五列本发明配准结果括号内的参数依次为：待配准图像H和最终配准图像H″之间实际横坐标平移量、实际纵坐标平移量、实际旋转角度。
[0070] 表2为采用粒子群(PSO)方法、鲍威尔(Powell)方法、仅采用点火次数矩阵M的配准方法以及本发明的配准时间的对比。从表中可以看到仅采用点火次数矩阵M的配准方法以及本发明在配准结果中时间上明显少于PSO(粒子群)方法和鲍威尔(Powell)方法的配准时间。同时本发明的配准时间比仅采用点火次数矩阵M的配准方法的时间更少，从而进一步的加快了配准速度。
[0071] 表2
[0072]
附图说明
[0073] 图1(A)为参考图像；
[0074] 图1(B)为待配准图像；
[0075] 图1(C)为配准后图像与参考图像间的差值；
[0076] 图2为本发明的流程框图；
[0077] 图3为计算边缘像素响应矩阵Y[n]各元素值Yij[n]基本公式；
[0078] 图4(A)笛卡尔坐标平面中参考图像的1～40个重心特征点；
[0079] 图4(B)笛卡尔坐标平面中待配准图像的1～40个重心特征点；
具体实施方式
[0080] 以下结合附图和实施例对本发明进一步说明。
[0081] 本发明实施例采用的是医学MRI(核磁共振成像)图像，如图1(A)、图1(B)所示，图1(A)为参考图像，图1(B)为待配准图像；。
[0082] 如图2所示，本发明实施例包括以下步骤：
[0083] 一.边缘检测步骤：
[0084] 由于MRI图像边缘清晰，采用sobel算子已能很好的寻找出图像目标的边缘，因而本发明实例采用sobel算子进行边缘检测。
[0085] 二.提取特征向量步骤：对图像进行ICM处理，得到特征矩阵：
[0086] 本实施例中，运行次数N为40，状态衰减系数f＝0.2，阈值衰减系数g＝0.9，阈值常数h＝20，分别对边缘像素采用原始ICM公式计算像素状态值，对非边缘像素采用改进的ICM公式计算像素状态值，图3为计算边缘像素响应矩阵Y[n]各元素值Yij[n]基本公式；
[0087] 得到点火次数矩阵Mij[n]，依次运行40次，计算每次运行得到的点火次数矩阵的重心得到重心特征点，如图4(A)、图4(B)所示，图中序号依次为每次运行得到的重心特征点的序号，最后得到40×2的待配准图像重心特征矩阵 和参考图像重心特征矩阵[0088] 三.计算配准参数步骤：
[0089] 计算H和R之间的x坐标平移参数Δx、y坐标平移参数Δy以及旋转参数Δθ；
[0090] 四.配准步骤：依据配准参数对待配准图像进行空间变换，实现图像配准；
[0091] 配准结果如图1(C)所示，图1(C)中的每个像素的灰度值越小，表示配准的结果就越好。
