混合样条大尺度医学图像弹性配准算法研究

混合样条大尺度医学图像弹性配准算法研究作者：天天论文网日期：2016-1-59:04:22点击：0摘要：针对大尺度形变医学图像配准速度慢和精度低的特点，提出一种结合薄板样条(TPS)和B样条的弹性配准方法。该方法采用尺度不变特征变换算法(SIFT)进行图像特征提取与匹配，利用TPS算法将特征点对作为输入进行预处理，以降低浮动图像的形变尺度，从而提高下一步B样条配准的速度与精度。然后使用局部区域细化层次B样条方法将TPS生成的较稀疏的形变网格作为初始网格，结合有限记忆优化算法(L-BFGS)对控制网格做进一步地处理，此过程只对形变较大的局部区域进行细化，以实现与参考图像的快速精确配准。实验结果表明，该方法较层次B样条方法有效地提高了配准的速度和精度。关键词：医学图像；弹性配准；薄板样条；层次B样条；大尺度形变引言：医学图像配准是医学图像融合、图像与标准图谱的匹配、外科手术规划与设计等研究的基础，对提高临床诊断治疗、病情监测、外科手术水平等有着非常重要的积极作用[1,2,3]。根据不同的临床需求，在配准时可采用不同的算法。在某些情况下，使用快速的刚性配准算法就能够达到要求，但是在浮动图像与参考图像之间存在大尺度形变的情况下，就需要通过非线性变换模型来描述复杂的形变过程从而解决问题，如形变复杂多样的软组织图像的配准、不同模态图像之间的配准、不同个体图像之间的配准等。针对这类问题，可将目前的方法大体分为基于径向基函数和基于伪物理模型两大类。前者采用径向基函数来表示形变域，如多项式函数[4]、薄板样条(ThinPlateSpline，TPS)[5,6]、B样条[7,8,9]、小波基函数等；后者使用偏微分方程来描述形变的平衡状态，如弹性模型、粘性流体模型、光流场模型等[10]。层次B样条[11,12,13]是目前常采用的弹性配准方法，它综合考虑了密集和稀疏的控制网格特点，力求在保证配准精度的同时，提高配准的速度。但其不足在于对每一层都要进行全局细化，不仅造成了大量不必要的计算，而且对于在上一层已经配准好的区域也将会受到整个形变域的影响。TPS算法[5]由于具有全局变形且没有冗长迭代的特点，同样是一种建模变形的有效工具。然而，TPS虽然可以实现快速计算，但是与层次B样条相比会产生较粗糙的配准结果。针对层次B样条配准算法很难对存在大尺度形变的浮动图像快速地达到理想配准精度的不足，本文结合TPS的快速和B样条的精度，采用TPS算法做预处理，然后使用局部区域细化层次B样条方法做进一步精确处理并与参考图像配准。实验结果表明，在配准的时间和质量上，该方法较层次B样条配准算法都得到了提高。1算法框架本文将TPS算法与局部区域细化层次B样条方法相结合，进行大尺度医学图像弹性配准。算法框架如图1所示。用函数)(XR和)(YF分别表示两幅输入图像，其中R为参考图像，F为浮动图像，X和Y分别代表各自图像的定义域。)(tS为相似度函数，其中t为控制参数，用均方差MSD(MeanSquareDifference，MSD)作为相似性测度来衡量两幅图像的匹配效果。此时，配准就是寻找一个形变函数，以使)(tS取得最小值。RR((XX))MMSSDDBB样样条条变变换换g(X)S(t)S*F(g(X))SIFT算法进行特征提取与匹配SIFT算法进行特征提取与匹配特征点对RR((XX))FF((YY))初步形变网格TTPPSS变变换换f(X)t计算所有控制点领域的MSD计算所有控制点领域的MSD插插值值提高局部区域的控制网格分辨率提高局部区域的插值控制网格分辨率插值LL--BBFFGGSS图1算法框架首先，输入待配准图像)(XR和)(YF，使用尺度不变特征变换算法(ScaleInvariantFeatureTransform，SIFT)[14,15]进行图像特征提取与匹配；将上一步得到的特征点对坐标作为TPS算法的输入，对X进行TPS变换，从而得到新的区域)(Xf坐标，通过插值生成初步形变网格；将该形变网格作为局部区域细化层次B样条的初始网格，对X进行B样条变换，从而得到区域)(Xg坐标，通过插值得到浮动图像在区域)(Xg的取值))((XgF；利用MSD计算参考图像)(XR和插值图像))((XgF间的相似度)(tS；将)(tS输入到优化模块中，利用有限记忆优化算法(LimitedmemoryBroyden-Fletcher-Goldfarb-Shanno，L-BFGS)[16]进行计算得到最优变换参数t，此过程需要通过迭代来实现，直到在该层使)(tS取得最小值；然后，对形变较大的局部区域，提高控制网格分辨率进行再一次配准，直到满足配准要求；最后，输出配准结果图像*S。2TPS算法降低形变尺度为了能够同时提高配准的速度和精度，本文采用TPS算法做预处理，而SIFT算法提取图像特征点并匹配是进行TPS预处理的前提。一般情况下，SIFT算法会产生大量的、健壮的匹配点对，根据经验，总共选用200对左右的匹配点对，就能够使TPS算法获得较好的初步形变网格。TPS算法是将插值问题模拟为一个薄金属板在点约束下的弯曲变形。它首先确保配准点集的一一对应，然后利用最小化TPS弯曲能来联合求解该点集之间的匹配矩阵及映射参数，使图像能够在控制点上达到精确吻合，并且在其约束下通过插值以使其他点也尽可能获得较好的纠正。表示浮动图像变形能量的代数式可定义为：22z(x,y)U(r)rlogr(1)它的解为UUxy22()2222，其中22yxr，)(rU是构建薄板样条的基函数。通过SIFT算法得到的浮动图像上的特征点坐标记为(x,y)。要使金属板在点(xi,yi)处的高度为iz，并且让其具有最小弯曲能量，即在二维空间中找到一个形变函数),(yxf，使薄板样条插值的罚函数)(fEs最小。该罚函数可定义为：Efdxdyyfxyfxfs()(()2()())22222222(2)令),(iiiyxP，ni,,2,1，设),(yxPri，为控制点对之间的欧几里德距离。定义3n矩阵1122111nnxyxyxyP(3)及nn矩阵121212120()()()0()()()0nnnnUrUrUrUrUrUrK(4)TKPL=P0(5)其中，0为3行3列的零矩阵，TP代表P的转置。构建行向量111222((,),(,),,(,))nnnVzxyzxy，z列x向y量(,0,0,0)TYV。通过下式：11(|)TxyaaaLYW(6)可以得到12(,,,)Tn。进而构造函数f(x,y)：nifxyaaxxayywiUPixy1(,)1((,))(7)该函数对于所有i有f(xi,yi)zi，并使罚函数Es(f)最小。通过f(x,y)得到的形变网格，将作为局部区域细化层次B样条的初始网格。3B样条算法细化形变网格并配准3.1B样条初步细化采用B样条自由形变模型[12]将TPS算法生成的较稀疏的形变网格覆盖在浮动图像上，通过调整控制点的位置，带动所嵌套的浮动图像发生形变。B样条具有良好的局部控制特性，浮动图像上某点处的位移值仅与其相邻的24个控制点有关。因此，图像中部分控制点移动仅会影响其周围有限范围内的空间变化，而不会影响整幅图像的像素点位置。__由于三次B样条具有二阶连续性和计算高效性等特性，因此本文采用三次均匀B样条基函数。为了能够根据TPS算法得到的网格控制点的位移值求得任意位置处的值，在二维空间中构造一个形变函数),(yxg，将与控制点相关的信息尽可能光滑地传播于浮动图像中的所有点。该形变函数的表达式为：mplqmlmgxyBluBv,3030(,)()()(8)其中1xp，1yqxxu，yyv)3,,0(lBl代表三次均匀B样条的l次基函数，其形式为：()(1)/63B0uu，()(364)/632B1uuu，()(3331)/632B2uuuu，()/633uuB，其中10u。它们代表控制点对),(yxg的权重值，即根据控制点到),(yx的距离来加权各个控制点对),(yxg的贡献。形变函数),(yxg生成的曲面分别由x方向上的3r个节点和y方向上的3s个节点同时控制。记TPS，每个控制点间的距离为ji,，它控制着形变程度。将浮动图像中点的运动位移分解为x和y两个方向，并分别求出所有点在这两个方向上的运动位移。各控制点在x和y方向上的位移值即为待优化的参数，网格中的每个节点有2个参数表示这个节点的位移，因此参数总数为2(r3)(s3)，通过最优化过程求得。网格中的每一个节点都需要迭代计算，迭代时间取决于L-BFGS优化算法和公差要求。通常要找到最佳的变形参数设置需要较长的计算时间。与TPS相比它虽然花费了更长的时间，但得到了一个更精确的结果。3.2局部区域精确细化当浮动图像的局部存在较大形变时，如果直接采用较小的控制网格进行配准，很可能会导致发生局部扭曲的现象。其原因在于需要优化的参数多且变动范围大，寻优结果很容易陷入局部极值。层次B样条采用的是由粗到细的控制网格，每一层的控制点间距由上一层控制点的间距二分而得，并利用上一层的配准结果进行更精确的形变配准，可以较好地解决这类问题。但是，层次B样条的每一层控制网格都是针对整幅浮动图像的，即每层都要进行全局细化。当浮动图像在初步配准后，大部分区域已经达到了较好的配准结果，只有局部区域还存在一些差异时，层次B样条方法不仅会造成不必要的计算，且对于已经配准好的区域也将会受到整个形变域的影响。因此，本文提出局部区域细化层次B样条方法来解决上述问题。在使用TPS生成的较稀疏的控制网格配准后，对两幅图像中所有控制点所对应的像素点的8领域计算MSD，对大于给定阈值的控制点领域，锁定该控制点所在的（4343）控制网格，并对此区域内其他点不再进行计算。当与该区域相邻的控制网格也含有领域的MSD大于所给阈值的控制点时，将控制网格合并。然后，提高这些需要进行进一步细化的局部区域的控制网格分辨率，以完成精确配准，如图2所示。在所有局部区域中使用同样的方法，直到所有的控制点的领域都满足要求时，结束配准。.图2局部区域精确细化4算法流程本文针对大尺度医学图像弹性配准问题，提出了以TPS算法做预处理，降低浮动图像的形变尺度，然后使用局部区域细化层次B样条方法进行精确配准的算法，解决了大尺度形变图像配准速度慢且精度低的问题。算法流程如图3所示：开开始始SSIIFFTT算算法法进进行行图图像像特特征征提提取取与与匹匹配配结结果果是是否否满满足足要要求求结结束束YNTTPPSS算算法法得得到到初初步步形形变变网网格格BB样样条条算算法法对对形形变变网网格格做做进进一一步步形形变变计计算算所所有有控控制制点点领领域域的的MMSSDD误误差差是是否否大大于于给给定定阈阈值值YN得得到到配配准准结结果果图图像像提高局部区域控制网格分辨率提高局部区域控制网格分辨率图3算法流程图配准的基本步骤如下：1)用SIFT算法提取R、F的特征点并匹配，得到特征点对坐标；2)将1)中得到的特征点对坐标作为TPS算法的输入，根据式（7）求得形变函数f(x,y)；3)将2)作为局部区域细化层次B样条的初始网格，根据式（8）结合L-BFGS优化算法计算变换参数；4)使用相似性测度函数判断参数是否达到最优，若是，则得到在较稀疏的控制网格下的最优解，否则，返回3)，计算新的变换参数，直到在该层使MSD取得最小值；5)计算两幅图像中所有控制点的领域的MSD，对于大于给定阈值的领域，提高该控制点所在的（4343）控制网格的分辨率，返回3)，计算新的变换参数，直到在该层使MSD取得最小值。在所有局部区域中使用同样的方法直到所有的控制点领域都满足要求时，结束配准。5实验结果及分析5.1图像配准结果在本实验中，使用本文提出的配准方法对大量的医学图像进行了仿真实验验证，图像大小均为512512，选取控制点领域的MSD阈值为80。实验平台为Intel(R)Core(TM)i3-2120CPU(双核四线程，3.30GHz)，4G内存。在MatlabV8.0环境下编写应用程序。图4和图5分别显示了两组肺部