一种基于椭球模型的非刚性医学图像抗旋转配准方法与流程

本发明属于图像处理技术领域，特别是医学图像配准技术领域，具体涉及一种基于椭球模型的非刚性医学图像抗旋转配准方法。

背景技术：

图像配准(imageregistration)就是将不同时间、不同传感器(成像设备)或不同条件下(天候、照度、摄像位置和角度等)获取的两幅或多幅图像进行匹配、叠加的过程，它已经被广泛地应用于遥感数据分析、计算机视觉、医学图像处理等领域。

图像配准是图像分析中的一项非常重要的技术，主要通过寻找某种变换，使两幅图像的对应点达到空间位置的一致，在医学诊断过程中，由于存在不同模式图像表现不同性质的物理机制、患者摆位的差异、成像参数的变化、不同成像设备间空间分辨率不同等现实问题，所以单单凭借医生手动将两张或两组不同模式的图像在空间上做对准受到很多局限，且常带有较大的主观性，其可靠性往往不高，不可避免地会产生误差。特别是在定向放射外科和颅脑手术可视化等应用领域，对图像配准的精度要求很高，使得医学图像配准成为一项必要而又相当困难的任务。

配准技术的流程如下：首先对两幅图像进行特征提取得到特征点；通过进行相似性度量找到匹配的特征点对；然后通过匹配的特征点对得到图像空间坐标变换参数；最后由坐标变换参数进行图像配准。

图像配准技术的基本过程主要分为三个步骤:

①寻找图像中的对应特征量并提取出来；

②根据特征量寻找最佳匹配变换；

③利用得到的最佳变换对待配准图像进行变换和参考图像进行匹配。

其中第二步是配准过程的关键，也是配准算法研究的核心内容。

当前常用的基于b样条的ffd模型,可应用于乳腺增强型mri图像的非刚性配准,但此方法因为采用互信息的配准度量并不适用于mri外某些图像，同时计算耗费时间相对较长；流体映射模型,用两组反映待配准图像轮廓信息的对应控制点集,对头颅,肝,脾等各种器官进行非刚性匹配，但如果不能实现流体模型的控制点的严格一致,图像配准有一定的失真；基于“de-mons-base”的算法,在判断出待配准图像上各个像素点运动的基础上，通过对各个象素点的移动来实现基于非刚性配准，但这种方法仅能在小参量下实现非刚性配准，当待配准图像有一定的旋转角度时，配准效果就比较差了。

如何建立合理的变形模型，适应各种组织变形，尤其是如何解决图像的旋转和平移问题，是图像处理研究的热点和前沿领域，本发明在总结前期非刚性领域已有的配准、融合研究工作经验的基础上，提出一种全新的改进模型——椭球模型，改变传统近刚性的组织图像的配准方法，解决特征点自动选取后、图像大旋转问题，提高现有算法鲁棒性差的问题，提高融合精度，并实现多模医学图像配准、融合与远程医疗信息系统的结合。

技术实现要素：

为了克服现有的医学图像配准方法存在的抗旋转差，不适应大角度变换，和配准速度慢等问题，本发明提供一种基于椭球模型的非刚性医学图像抗旋转配准方法。所述方法在获取图像特征点和最小特征矩形域的前提下，采用椭球模型对图像的最小特征矩形域进行划分，求取每个划分内像素灰度的平均值，作为该区域像素灰度的代表；依托像素灰度序列，采用傅里叶变换算法，求取图像旋转变换参数；对图像进行一次旋转变换，如此完成一次配准过程；然后递减图像划分的子区域块数，重复上述过程，直至图像划分块数为1，此时对图像进行插值运算，实现图像配准。所述方法与一般的配准方法比较，减少了搜索位置，大幅提高了配准速度；抗旋转性好，配准更为准确。

为实现上述目标，本发明采用以下技术方案：

一种基于椭球模型的非刚性医学图像抗旋转配准方法，所述方法在对图像进行初始变换，获取图像特征点和最小特征矩形域的前提下，采用椭球模型对图像的最小特征矩形域进行划分，求取每个划分内像素灰度的平均值，作为该区域像素灰度的代表；依托像素灰度序列，采用傅里叶变换算法，求取图像旋转变换参数；对图像进行旋转变换，完成一次配准过程；通过减少椭球模型的划分参数，对图像进行再次划分和变换，逐步提升图像配准精度；重复上述过程，直至图像划分块数为1，此时对图像进行插值运算，实现图像配准。

一种基于椭球模型的非刚性医学图像抗旋转配准方法，所述方法包括以下步骤：

1)、图像基本变换：对原始图像进行某种变换，其目的就是为了求取图像的特征点；所述基本变换方法包括但不限于“基于小波变换的边缘检测”，所述“基于小波变换的边缘检测”就是将原始图像进行小波变换，将其分解在不同频段，找出高频部分模的极大值后，进行筛选从而得到图像特征点。

2)、构建图像的包含最有价值特征点的最小矩形域：构建最小矩形包含图像的特征点，以此作为确定待配准图像与参考图像的包含最有价值特征点的最小矩形域。

3)、以椭球模型对图像的包含最有价值特征点的最小矩形域进行划分：以特征矩形域模板的中心为对称中心，以a(n-i)/n(i＝0,1,2,...,n-1)和b(n-i)/n(i＝0,1,2,...,n-1)为长短焦半径，画同对称中心等离心率的椭圆，这样在特征矩形域的基础上，画出了n个椭圆环，分别记作ci(i＝0,1,2,...,n-1)。

以特征矩形域模板的中心为圆心，并作为起始点，设定起始位置，对每个椭圆环再进行划分，将椭圆环ci按角度均分成n个扇形区域，这样就将椭圆形模板划分为n²个扇形区，并将椭圆环内每个扇形区按顺时针方向做好标记，将其记为n²个扇形区sij(i,j＝0，1，2,...,n-1)。

将每个扇形区内截取的环作为一个单位，每个单位对应的灰度值为该扇形区域内所包含的像素灰度的平均值fij。这样，每一圆环便输出n个灰度值，并将其构成一维数组fn(x)。靠近圆心的圆环包含像素少，中心部分用扇形区代替扇形环。

4)、确定旋转系数，实现初次旋转：对旋转前和旋转后得到的一维数组分别进行傅立叶变换，求其幅度谱，根据傅立叶变换的平移不变性，可知一维数组平移前和平移后其幅度谱是不变的，作为图像匹配的旋转不变量，用傅立叶变换的幅度谱进行相关计算，寻求各相关系数的平均值的最大值，这就是确定第一组旋转系数fn，对待配准图像进行第一次旋转变换得到图像in。

5)、多次重复旋转：减小椭球划分模型的参数n，例如但不限于取n＝n/2，对in待配准图像进行椭球模型拆分，重复步骤3)和4)，确定第二组旋转系数fn/2，对待配准图像in进行旋转变换得到图像in/2。重复步骤5)，直至n＝1。

6)、当n＝1时，对待配准图像进行最后一次旋转变换得到图像i1，对图像进行插值等处理后即得到配准后图像。

所述步骤2)中构建包含最有价值特征点的最小矩形域方法包括以下步骤：

a)、对图像初始变换如小波变换中确定的特征点进行判别，如果确是特征点，则将其加入特征点集合a中，构建待选点集合b，b初始为空集；

b)、选取特征点a1和a2作为左上角和右下角，构建矩形rn-1，将a中选中的点移入b中；

c)、判断集合a是否为空集，如果集合a为空集，则得到最小矩形域，算法结束；如果集合a不为空集，则进入下一步；

d)、从集合a中选择点ai并将其移入集合b中；

e)、选择特征点ai和集合b中的另一点作为左上角和右下角，构建矩形rn-1，返回步骤c)。

所述计算图像旋转参数的方法包括以下步骤：

a)、将最小矩形区域按照椭球模型拆分成同中心椭圆环；

b)、按角度均等拆分椭圆环为扇形区域；

c)、计算每个扇形区域的像素灰度平均值fij；

d)、通过傅里叶变换计算，求取各相关系数平均值的最大值fn。

本发明的优点和有益效果为：

1)、本发明所述方法与一般的配准方法比较，减少了搜索位置，从而大大提高了配准速度，对于一幅n×m个像素的图像，采用全像素配准算法时，理论上需要的搜索时间为o(nm)次，而本发明所述算法需要的搜索时间为o(lognm)次(o(x)是随x的值单调变化的时间度量)，扣除预处理构建扇形区域图像所花费的时间，整个配准过程所耗费的时间仍减少了很多。

2)、本发明所述方法与普通旋转配准方法比较，不需要特征点的完全一致对应，只需要构建出最小特征矩形域即可实现配准，应用更广泛；在任意角度旋转时进行匹配，进行加权修正后的算法对基准图和模板图之间存在灰度、噪声、光照和对比度等变化有较好的适应性，配准更为准确，抗旋转性好。

附图说明

附图1是本发明所述方法的流程图。

附图2是本发明所述的最小特征矩形域构建方法流程图。

附图3是本发明所述的通过线性函数计算图像旋转角度的流程图。

附图4是本发明所述方法实施方式示意图。

附图5是本发明所述的按照椭球模型对特征矩形域进行首次划分的示意图。

附图6是本发明所述的按照椭球模型对特征矩形域进行再次划分的示意图。

附图7是本发明所述的按照椭球模型对特征矩形域进行倒数第二次划分(n＝4)的示意图。

附图8是本发明所述实施例中的原始参考图像。

附图9是本发明所述实施例中的原始待配准图像。

附图10是本发明所述实施例中选取n＝64时构建的原始参考图像的椭球模型图像。

附图11是本发明所述实施例中选取n＝64时构建的原始待配准图像的椭球模型图像。

附图12是本发明所述实施例中选取n＝64时经傅里叶旋转变换后的配准图像。

附图13是本发明所述实施例中选取n＝16时构建的原始参考图像的椭球模型图像。

附图14是本发明所述实施例中选取n＝16时构建的原始待配准图像的椭球模型图像。

附图15是本发明所述实施例中选取n＝16时经傅里叶旋转变换后的配准图像。

附图16是本发明所述实施例中选取n＝4时构建的原始参考图像的椭球模型图像。

附图17是本发明所述实施例中选取n＝4时构建的原始待配准图像的椭球模型图像。

附图18是本发明所述实施例中选取n＝4时经傅里叶旋转变换后的配准图像。

附图19是本发明所述实施例中选取n＝1时经傅里叶旋转变换后的配准图像。

附图20是本发明所述实施例中选取n＝1时经简单融合后的效果图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

参见附图4、附图5、附图6、附图7，本发明的实施方式是，首先对参考图像和待配准图像分别进行变换，提取图像特征点，进而得到参考图像和待配准图像包含最有价值特征点的最小特征矩形区域；其次，运用本发明所述的椭球模型对参考图像和待配准图像进行拆分，将其拆分成扇形区域；

第三，运用本发明所述的傅里叶变换方法，计算图像旋转变换的参数；

第四，对待配准图像进行旋转操作；

第五，减少椭球模型中的拆分参数，重复步骤二至步骤四，直至得到最大拆分区域，即n＝1；

第六，对待配准图像进行插值处理，完成图像配准工作。

实施例

参见附图8至附图20。选取visualhuman图像库中的标准图像ct1109作为原始参考图像(附图8)，图像pd1110作为原始待配准图像(附图9)。

对图像ct1109和pd1110分别进行基于小波变换的边缘检测操作，提取图像的特征矩形域后，选取n＝64，以椭球模型构建的基准图像和待配准图像分别如附图10和附图11所示，经傅里叶旋转变换后的配准图像如附图12所示。

选取n＝16，以椭球模型构建的基准图像和待配准图像分别如附图13和附图14所示，经傅里叶旋转变换后的配准图像如附图15所示。

选取n＝4，以椭球模型构建的基准图像和待配准图像分别如附图16和附图17所示，经傅里叶旋转变换后的配准图像如附图18所示。

选取n＝1时，以椭球模型构建的经傅里叶旋转变换后的配准图像如附图19所示，经简单融合处理后的配准图像如附图20所示，由附图20可知，通过本发明所述的方法，实现了图像的精确配准。

最后应说明的是：显然，上述实施例仅仅是为清楚地说明本发明所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本发明的保护范围之中。