医学序列图像（精选7篇）
医学序列图像 篇1
摘要：针对医学序列图像的特点,概述了两种压缩的方法[1-3],并在此基础上提出了一种基于最大互信息配准法的压缩方法,更有效的减少序列图像之间的相关冗余。该方法充分利用了帧间相关性,具有很好的压缩性能,可以有效满足现代医学图像的压缩要求。
关键词：医学序列图像,最大互信息配准法,研究
1 图像压缩。
所谓医学序列图像就是对某一相同的部位从相同的角度进行连续拍摄所得到的许多图像,其中的每一幅图像称为一帧,通常的医学序列图像都具有十几帧或二十多帧,甚至更多。这些图像的不同之处,主要在于病变区,而非病变区的信息基本相同。还有一种准序列图像,它是在一段时间内,对相同部位,从相同角度进行拍摄或者是从不同角度对同一部位进行连续拍摄,虽然角度不同,但拍摄的部位有相当一部分是重合的。从图像的整体来看它们也具有许多相似点,我们可以充分利用这些相同和相似处,找出图像之间的冗余信息。针对图像各帧之间存在很大相关性,采用了基于最大互信息法配准的技术,即寻找到一种变换,使得两帧待配准的图像间的互信息最大,从而尽可能地减少冗余信息,提高压缩效果。
2 压缩方法
2.1 差值序列压缩
所谓差值序列压缩,其实质就是一种预测方法,用前一帧作为下一帧的预测。就是除了第一帧要进行帧内压缩外,压缩其它的帧就是压缩该帧与其前一帧的差。因为在医学序列图像中,除了病变区,周围的图像及其象素值基本不变,如果对相邻两帧求差值,那么将能够得到不少的0值,这对压缩非常有利,整幅图像的压缩比也会有较大提高。
2.2 基帧选择压缩
上述序列图像的差值压缩,是对相邻的两帧求差值,然后压缩;压缩顺序是从第一帧开始,接着压缩第二帧,依次下去,直到压缩最后一帧,第一帧采用无损压缩,压缩第二帧就是压缩第二帧与第一帧的差值,压缩第M帧就是压缩第M(M>1)帧与第(M-1)帧的差值。但相邻两帧的差值压缩效果并不一定最好,可能帧M与不相邻的某帧N的差值压缩效果会更好。因此,压缩帧M时可以在已压缩的帧中选择一帧作为基帧,然后对它们进行差值压缩。对于差值图像,要提高其无损压缩的效率,有两点比较重要:
第一,差值图像中0的数目;
第二,象素值的位数的多少;如果象素值的位数都比较少,而0的数目又比较多的话,对提高压缩效果很有帮助。
由于以上两点,这里的基帧选择方法采用绝对误差准则:即求两帧的象素差值绝对值的平均值MAD,与将被压缩帧的绝对误差最小的已被压缩帧,将被选为基帧。
X(i,j)表示第x帧的第i行第j列的象素值。
基本算法如下:
a.用压缩算法压缩第一帧;
b.压缩第二帧(压缩第二帧与第一帧的差值),M=2;
c.压缩第M+1帧(M>=2):用绝对误差准则从已经压缩的1-M帧中选择一帧作为基帧,然后将第M+1帧与该基帧求差值,进行压缩。直到压缩完最后一帧。采用选择基帧的方法,能使压缩的效率进一步提高。
重点在于选择了基帧之后,采用最大互信息法配准技术[[44]],就是在基帧和待压缩帧之间寻找到一种变换,使得两帧待配准的图像间的互信息最大,从而尽可能地减少冗余信息,提高压缩效果。
3 基于最大互信息法配准两帧图像
3.1 互信息的计算
互信息是信息理论中的一个基本概念,通常用于描述两个系统间的统计相关性,或者是一个系统中所包含的另一个系统中信息的多少。对于待配准的两幅源图像A(x,y)和B(x,y),其互信息(Mutual Information,MI)定义如下:
式中,L为图像总的灰度级数,PA,B(i,j),为图像A、B中像素灰度的联合概率分布,PA(i)和PB(j)分别为图像A和B中像素灰度的边缘概率分布。PA,B(i,j),PA(i)和PB(j)可以由联合直方图函数hA,B(i,j)归一化求得:
联合直方图hA,B(i,j)是通过统计两幅图像的对应位置的灰度对(i,j)出现次数得到的。
3.2 配准算法
将互信息用于图像配准的关键思想是,如果两幅图像达到匹配,则它们的互信息达到极大值。整个配准过程如图1所示。配准过程为:首先对配准图像进行必要的预处理,参考图像始终保持不变,对浮动图像进行图像变换,与参考图像比较,通过联合直方图统计得到对应的互信息,通过最优化搜索算法,搜索得到是互信息达到最大时对应的配准参数,根据这一组配准参数对浮动图像进行变换即可与参考图像达到匹配。
3.3 优化算法
把互信息MI作为测度函数,当互信息MI达到最大值时,就认为已经找到了最佳参数。在优化方法中,通常是对代价函数求最小值,这时采用将MI的公式前加上一个负号,即可作为代价函数。基于最大互信息的配准过程实际上是一个寻找多参数的优化过程。在中采用了刚性配准,进行沿X轴方向的偏移量,沿Y方向的偏移量及沿图中心的旋转角度三个参数的最优化搜索。搜索方法采用POWELL算法,它是在研究具有对称正定的二次函数(公式5)极小化问题时形成。
其基本思想是在不用导数的情况下,在迭代法中采用精确一维搜索的条件下逐次构造Q共轭方向。因此,Powell法本质上属于共轭方向法。而共轭方向具有二次收敛性,因此Powell方向搜索法有较高的收敛速度。关键步骤:求出沿X轴方向的偏移量,沿Y方向的偏移量及沿图像中心的旋转角度三个参数,经过三轮迭代后,得到了三个Q共轭向量P1、P2、P3。对于只有三个变量的正定二次函数,经过3次直线搜索必收敛于正定二次函数,即可求得极小点。
4 结论
中的序列图像压缩算法,均可达到一种无损压缩,同时能提高医学序列图像的压缩比。虽然,在现实中,医学序列图像的非病变区并不是完全一样,但其仍然比较相似,对无损压缩仍能起到相当的帮助作用。
参考文献
[1]周志勇.JPEG2000静止图像压缩算法研究[D].四川大学硕士学位论文,2004,5.
[2洪英杰.新一代静止图像压缩算法JPEG2000的研究与实现[D].南京理工大学硕士学位论文,2004,6.
[3]陈静.基于JPEG2000的医学图像压缩研究[D].泰山医学院硕士毕业论文,2004,5.
[4]高智勇,顾滨.林家瑞.基于互信息的医学图像配准实验.生物医学工程学杂志,2003,20(3):476-478.
[5]傅英定,成孝予.最优化理论与方法[M].北京:电子科技大学出版社,1996.
基于序列图像的三维重构 篇2
利用一系列已经标定好的二维图像精确的恢复物体的三位模型，是目前机器视觉领域研究的一个热点。其广泛应用于遥感图像分析、机器人视觉、虚拟现实、精密工业测量及军事侦察等领域。通常的双目立体视觉是建立在两幅图像上的重构，局限在于不能得到完整的三维模型。但通过近十多年的发展，双目立体视觉重构技术已经十分成熟。
本文就是在相机标定的基础上[1]，利用对极几何约束的方法，先将两幅图像按照双目立体视觉进行重构。然后利用剩余图像与重构模型间的匹配关系，采用反投影映射的方法计算出剩余图像的外方位元素。通过较高精度的相机标定及重构中对各参数的迭代优化控制累计误差。
1 基于对极几何重构
对极几何关系是指从两个不同视点获取的来自同一个场景的两幅图像之间的约束关系，这种约束关系可以用基本矩阵(Fundamental Matrix)来描述。基本矩阵包含了相机从第一个位置到第二个位置相对运动的所有信息。
1.1 相机模型
机器视觉中常用的成像模型是透视相机模型。满足投影方程：
其中，m为二维图像坐标m=[x y1]T，相应的，M是空间中三维点的坐标M=[X Y Z1]T,P是3×4的投影矩阵。在欧式重构中P可以表示为：
其中：
纵横比，s表示畸变因子，[u v]为图像主点也就是相机中心在图像中的坐标。通常s=0，α=1，图像主点也接近于图像的中心。R和t为相对于世界坐标系的旋转矩阵及平移矩阵。
1.2 两幅图像的重构
基本矩阵F是3×3秩为2的矩阵，具有7个自由度。该矩阵包含了两幅图像间的变换关系，因而对基本矩阵F的稳健估计是三维重构的基础。如图1所示，三维空间中某个点M在图像1和图像2上的投影分别为点m和m'，这两个点即为一组匹配点。基本矩阵与匹配点满足下面关系式：
由式(3)可知，基本矩阵可以采用RANSAC的方法由两幅图像的匹配点求得[2][3]。
由对极几何关系[4]可知：F矩阵与投影矩阵P有如下关系：
由PC=0求得。假设以第一幅图像的主点为世界坐标原点，可得两幅图像的投影矩阵为：
将式(5)代入式(4)可得：
其中，E称为关键矩阵，可由内方位元素矩阵K和基本矩阵F求得。对式ET t=0奇异值分解可以得到t的单位向量，由方程RT ET=tT可得旋转矩阵R。
将两幅图像的投影矩阵和匹配点代入投影方程式(1)可以求得一个基于两幅图像的三维初始模型。
2 序列图像重构
在得到两幅图像重构结果的基础上，本文采用每次将一幅图像加入原有模型的方法完成整个模型的重构。首先，根据新加入图像与初始的三维模型间的匹配关系确定该图像的外方位元素，然后将新产生的三维坐标点加入原有模型[5]。
如图2所示，对于新加入的图像上的点mi，其和前一幅图像间的匹配点mi-1可由RANSAC方法求得，若点mi-1对应的三维空间点为Mi，则新加入图像上的点im所对应的空间点也为Mi。同理可求新加入图像上所有的匹配点所对应的三维点，将所有对应的二维点集[m1 m2 m3]和三维点集[M1 M2 M3]代入投影方程式(1)便可以得到新加入图像的投影矩阵P
由于误差和噪声的存在，每加入一幅图像，采用如下的优化函数对其投影矩阵进行校正：
3 稠密重构及纹理贴图
为了使重构结果的表面信息较为完整，本文采用Quasi稠密匹配的方法[6]求得致密的匹配点，同时经过序列图像的重构所得到的投影矩阵重构出大量的三维空间点，最后利用delaunay三角剖分法对重构结果进行分割及贴图。
4 实验结果及分析
仿真实验中，相机的内方位元素经标定为：
随机拍摄得到五幅大小为1024×768的序列图像,如图3所示。
先采用SIFT匹配的方法[7]求得初始匹配点，利用RANSAC的方法剔除误匹配的同时，得到每两幅图像之间的基本矩阵。选择匹配数目最多的两幅图像进行初始的重构。最后运用本文基于序列图像重构的方法对该序列图像重构，结果如下图所示。
由表1可以看出，由于第三幅图像和第二幅图像间的视角变换比较大，导致匹配点数目较少，基本矩阵的精度较低，从而使反投影误差较大。由于采用了式(8)的的函数进行优化，第四幅图像和第五幅图像的反投影误差并没有明显的增长。因此，在相机标定精度较高的情况下，本算法并不会由于中间图像出现较大误差而影响后续图像的重构。
5 结论
本文提出了一种基于序列图像的三维重构算法。在相机标定的基础上先用两幅图像进行初始重构，通过每次添加一幅图像的方法求得剩余图像的外方为元素，采用稠密匹配的方法得到大量的匹配点，结合所有图像的外方位元素重构出完整的三维模型。从实验结果来看，在相机标定精度较高的情况下，采用该算法能够得到较低失真的三维模型。
实验结果表明该算法仍有不足的地方：
1)该方法得到大量的三维空间点有一部分是重复的点，应该剔除;
2)由于稠密匹配耗费时间比较长，因此实验总耗时比较大，对稠密匹配需要进一步研究;
3)该算法受限于相机标定的结果，特别对于图像有畸变的情况，应先对图像进行校正。
摘要：提出了一种基于序列图像的三维重构方法。在相机标定的基础上,利用两幅图像间的匹配关系得到一个初始的三维模型,然后利用剩余图像与初始模型间的匹配关系,采用反投影映射的方法计算出剩余图像的外方位元素。最后经过稠密匹配、三角剖分和纹理渲染得到物体的三维模型。实验结果表明,相机标定精度较高的情况下,采用该算法能够得到较低失真的三维模型。
关键词：三维重构,对极几何,外方位元素
参考文献
[1]Zhang Zhengyou.A Flexible New Technique for Camera Calibation[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2000,22(11):1330-1334.
[2]M A Fischler,R C Bolles.Random sample consensus:a paradigm for model fitting with applications to image analysis and automated cartography[J].Communication of ACM,1981,24(6):381-395.
[3]R Hartley.In defence of the8-point algorithm[C].In Proceedings of the Fifth International Conference on Computer Vision.1995:1064-1070.
[4]R Hartley,A Zisserman.Multiple view geometry in computer Vision[M].Second.New York:Cambridge University Press,2003:241-250.
[5]M.Pollefeys,L.Van Gool.Visual modeling with a hand-held camera[J].International Journal of Computer Vision,2004,59(3),207-232.
[6]J.Kannala,S.S.Brandt.Quasi-dense wide baseline matching using match propagation[J].IEEE Conference on Computer Vision and Pattern Recognition,2007.
心脏超声序列图像轮廓的弹性配准 篇3
心脏的运动过程中包含着大量的信息。根据所得心脏超声序列图像进行诊断时, 需要观察心脏特定部位的收缩和舒张运动。本文应用配准技术[1]对超声序列图像特定部位进行逐帧配准, 将它作为心脏序列帧图像跟踪研究的一种基础技术。
由于心脏是一种变形体, 心脏的运动是由功能性的舒缩等变形运动和呼吸、血液涌动和结构间牵引等引起的平移旋转等非功能性运动构成的复合运动[2]。因此为了得到满意的结果, 必须使用非刚体配准。
弹性变换[3]是非刚体变换中的一种特殊的变换, 它允许局部调整, 并符合某种连续性或平滑性要求。本文应用薄板样条函数进行点约束算法, 对心脏超声序列图像进行弹性变换, 并使用了弹性近似方法来解决存在标记点定位误差情况下的弹性变换模拟问题。
本文中超声序列图像弹性配准方法主要包含4个步骤: (1) 在待配准图像中选取一定数量的标记特征点; (2) 建立起两幅图像的标记特征点之间的对应关系; (3) 利用一种插值方法求取图像之间的配准变换; (4) 将求得的变换作用于待配准图像, 实现图像的弹性配准。
2、标记特征点的选取
在超声序列图像配准过程中, 标记特征点的选取和精确对应至关重要。它可以通过手动选点, 也可以自动、半自动确定。
由于手动选取标记特征点的方法费时费力, 自动选取的方法会导致局部失真扩散到全局, 因此我们采用了半自动标记特征点的选取方法。首先, 我们在基准图像中选取较能反映运动变化的, 并且具有较强的特征性, 容易与周围邻近点区分的点, 例如, 灰度局部极大值、局部边缘、角点、不同结构交叉点等作为标记特征点。然后, 以标记特征点为中选取一定大小的图像块作为匹配块, 由于心脏各结构部位的相邻时刻的运动范围具有一定的大小。因此, 我们可以在待配准的图像中, 标记特征点坐标位置附近的某一范围内, 搜索得到相关性最大的点为标记特征点的对应点, 采用这种方法获取得到6组特征标记点对。
3、薄板样条平滑算法
无论是哪一种标记特征点的选取方法, 误差都是不可避免的, 因此考虑到对应点之间的误差, 采用一种平滑的处理方法, 对常用的弹性配准的插值方法薄板样条法进行平滑处理, 使得在配准时, 对标记特征点的精确性和数量的要求都明显降低。
3.1 薄板样条插值
薄板样条插值是一个多变量插值问题[4], 即从待配准的两幅图像中 (图象维数为d) , 标定n对对应标记点pi和qi, i=1, …, n, 在一个合适的Hilbert空间H上寻找连续变换u:Rd→Rd, 满足: (1) 最小化一个给定的泛函J (u) :H→R; (2) 实现如下插值
插值过程形象地模拟为一个薄金属板在点约束下的扭曲变形, 要使金属板在点 (xi, yi) 处的高度为zi, 并且该板具有最小弯曲能量, 即薄板函数u (x, y) 使泛函J (u) 最小。对于变换u的d个子变换u i, i=1, …, d, 泛函J (u) 可以分解为d个子泛函J (ui) 。因此求换u的问题可以分解为d个子问题, 即对于d维图象, 在每一维上分别求取ui。当d维和导数阶次为m阶时, 泛函可以表示为
其中单个泛函是
αk是正整数, x为坐标矢量, xk为x的单个分量。此泛函是旋转、平移不变的。当d=m=2时, 此泛函就是文献中应用的函数。泛函式 (公式3) 的零空间∏m-1 (Rd) , 是由Rd上的阶的数从0到m-1阶多项式函数Φν (ν=1, …M) 张成的空间。它的维数为:
M必须小于对应点的个数n, 这就决定了对应点的最小数量。d=m=2时, 对应点的数目必须大于3。最小化泛函式 (公式3) 的解析形式是:
U (x, pi) 为插值函数, 它是由维数d, 泛函的最高阶次m及Hilbert空间H决定的。选择Sobolev空间S=H2, 则
基函数U (x, pi) 张成一个n维函数空间, 此空间是仅由基准图像的标记点pi决定的。当d=m=2时, 得到二维图像的插值基函数U (, 其零空间是由函数Φ1 (x) =1, Φ2 (x) =x, Φ3 (x) =y张成。
解下面的线性方程组来求解方程式 (公式4) 的系数a= (a1, …, aM) T, 和w= (w1, …wn) T
Ki, j=U (pi, qi) , Pi, j=Φj (pi) , ν是待配准图像中的标记点qi在每一维坐标上的列向量, PWT=0是边界条件, 保证在无限远处弹性变换为0。
3.2 薄板样条的平滑
平滑薄板样条插值是在假定标记点的位置能准确选定的基础上进行计算的。而在实际中。标记点的位置只能近似选定, 存在误差。由于薄板样条插值对点位置的偏差敏感, 因此如果对应标记点的位置偏差较大, 就不能得到准确的配准结果。为此, 需对插值方法进行平滑处理, 以减少标记点偏差带来的影响。
引入一个二次项, 对泛函式 (公式2) 做如下处理
式中的第1项称之为数据项, 它是标记点变换后与目标点之间的欧氏距离之和;第2项为平滑项, 用于调节变换结果的平滑程度。最小化此泛函, 得到变换, 它能够: (1) 最大程度配准标记点集; (2) 配准变换足够平滑。变换的平滑程度是由参数λ决定的, λ为正数。当λ小的时候。得到一个平滑程度小的弹性变换;当λ大的时候, 得到一个平滑程度高而变形程度小的变换。因此, 要根据图象的实际情况, 选择合适的介于0～∞之间的λ值来实现理想的弹性配准。
而泛函式 (公式3) 和式 (公式7) 有着相同的基函数, 求解系数a= (a1, …aM) T, 和w= (w1, …, wn) T的方程组在形式上是相同的, 即只需将方程组式 (公式6) 中的矩阵K用K+λI代替
4、心脏左室短轴轮廓超声序列图像的配准
利用以上所述的方法, 我们采用半自动的方法选取标记点, 即在基准图像中人工选取6个点作为标记点, 并在待配准图像的相应位置的一定范围内, 通过相关性计算得到对应点, 并以这6组标记点采用薄板样条平滑算法进行插值, 对相邻帧的心脏短轴室壁图像进行配准, 并通过试验确定当λ=0.8。其结果如 (图一) 所示:
图 (a) 是在基准图像中提取的轮廓曲线。图 (b) 是将基准图像中的轮廓经过弹性配准变换后得到的轮廓曲线。图 (c) 是配准得到的相邻两帧图像的轮廓比较, 体现了相邻帧心脏短轴室壁的运动和变形, 说明了运用基于薄板样条平滑算法弹性配准方法进行轮廓跟踪的可行性。
5. 结束语
由于超声心动图所检测的心脏的运动, 是心脏各结构的功能性 (如舒缩变形等) 运动和呼吸、血液涌动和各结构运动的互相牵动等引起的非功能性运动 (如平移旋转运动等) 的复合运动。国际上所有的超声心动图运动检测如M型、组织多普勒及现有的全方向M型心动图[6]都存在这样一个问题:由于检测波束或方向线束是空间固定坐标, 为此取得的运动信息也只能是复合运动, 而不是只反映本结构机能的功能性运动, 所以这样得到的波形图不能准确反映心脏该结构某个部位的运动情况。因此我们尝试采用薄板样条平滑算法对超声序列图像进行配准, 取得了一定的效果, 并为心脏某结构部位的功能性运动与非功能性运动分离的提供了基础性研究。
摘要：心脏的运动过程中包含着大量的信息, 根据所得心脏超声序列图像进行诊断时, 需要获取心脏特定部位的收缩和舒张运动信息。本文尝试将基于薄板样条的弹性配准方法运用于心脏超声序列图像短轴室壁轮廓的配准, 从而获取其运动信息, 文中详细介绍了心脏超声序列图像的标记特征点选取和平滑薄板样条弹性配准的原理及应用结果, 该研究为心脏序列帧图像跟踪研究提供了一种新方法。
关键词：超声序列图像,薄板样条,弹性,配准
参考文献
[1].罗述谦, 吕维雪.医学图像配准技术.国外医学生物医学工程分册[M], 1999, 22 (1) :2-8.
[2].张莉, 吴剑锋, 林强.B超心脏序列图象左室短轴轮廓曲线的动态跟踪[J].福州大学学报 (自然科学版) .2000, 12.vol.28, No.6
[3].Toga A, Thompson P.The role of image registration in brain mapping.Image Vision Compute, 2001, 19:3-24
[4].Rohr K, Stiehl H S, Sprengel R et al.Point-based elastic registration of medical image data using approximating thin-plate splines[A].In:Proceed-ings Fourth International Conference on Visualizationin Biomedical Com-puter (VBC'96) [C].Berlin, Germany:Springer Verlag, 1996, 1131:297-306.
[5].罗述谦, 阎华.基于薄板样条的MRI图像与脑图谱的配准方法[J].中国生物医学工程学报, 2004, 23 (6) :479-485.
医学序列图像 篇4
手势生动、形象、直观,而且蕴涵着丰富的信息,是人与人沟通的重要媒介。手势跟踪则是在图像序列中确定出各帧中手势的位置、大小和方向等信息,是目前基于手势交互的重要内容。
目前,目标跟踪的算法[1,2,3,4]有很多,大体上可以分为基于匹配的跟踪和基于运动特性的跟踪。根据提取目标的特征不同,基于匹配的跟踪可分为区域匹配跟踪、轮廓匹配跟踪和特征匹配跟踪。基于运动特性的跟踪又分为基于光流的跟踪算法和运动预测跟踪算法。其中,Mean Shift算法因其计算简单、实时性好而得到了广泛应用。笔者将Mean Shift算法应用于序列图像中的手势跟踪,取得了较好的效果。
2 基于颜色直方图的Mean Shift算法
Mean Shift[5,6]是一种基于核函数密度估计的图像特征分析方法。在手势跟踪中,Mean Shift算法利用图像的颜色统计直方图作为整个搜索匹配的特征,并且采用梯度优化方法来减少特征搜索匹配的时间,实现快速的目标定位,同时利用Bhattacharyya距离作为对目标模板和候选目标的相似性函数,完成特征的匹配。
1)目标模板的表示
假设{xi}i=1,2,…,N表示目标模板区域中的像素位置,且目标区域的中心为x0,特征值bin的个数为m,则目标模板的特征值u=1,2,…,m的概率密度为
式中:k(x)为核函数,是一个各向同性的、凸的且单调递减的函数,它的作用是给目标模板区域的像素设置权值,即给远离目标模板中心的像素设置较小的权值,而给靠近目标模板中心的像素设置较大的权值。因为远离目标模板中心的像素容易受到其他目标或背景像素的影响,所以加权后的目标模板像素在进行密度估计时增加了稳健性,提高了搜索跟踪能力。函数k(x)中的的作用是为消除不同大小目标计算时的影响,将目标区域归一化。δ(x)是Delta函数,δ[b(xi-u)]总的作用是判断目标区域中像素xi的颜色值是否属于第u个bin,属于则值为1,否则为0。C是一个标准化的常量系数,使得。因此
2)候选模板的描述
设运动目标在第2帧及以后每帧中可能包含的目标区域为候选区域,设是候选区域的像素位置,当前帧中心坐标为y,候选模板的特征值u=1,2,…,m的概率密度为
式中:是归一化常数。
3)基于Bhattacharyya系数的相似性函数
引入Bhattacharyya系数[7]来度量候选目标和目标模板之间的相似性,其定义为
值在0~1之间,其值越大,表示两个模板越相似。在当前帧中不同的候选区域计算得到候选模型,使得最大的候选区域即为在本帧中的目标位置。
4)目标定位
为了在被检测帧中定位目标,应当最大。定位过程可以看作是从前一帧的目标中心y0开始,在其邻域内搜索目标的过程。将当前帧的目标中心定位为前一帧的目标中心位置y0,对式(4)在处进行泰勒展开,Bhattacharyya系数近似为
其中
式(5)中第1项是与y无关的,只有第2项随y变化,所以只需要对第2项进行分析。令
要使式(5)最大就是求式(7)最大,计算式(7)的Mean Shift向量,从候选区域中心y0递归地计算出新目标位置,即
式中:g(x)=-k′(x),Mean Shift算法是从起,向两个模型相比颜色最大的方向移动。
3 颜色概率模型
为了实现手势的跟踪,选取手势的颜色来作为跟踪特征,然后通过颜色概率模型把每帧图像转化为颜色概率分布图。常用的RGB色彩空间[8,9]并不能很好地反映视觉上的颜色感知性。由于笔者研究的是序列图像中的手势跟踪,因此这里的颜色就是肤色。由于肤色对HSV(色调、色饱和度和亮度)色彩空间的H分量很敏感,所以采用H分量的直方图建立肤色概率模型,即
式中:R,G,B分别代表RGB色彩空间中的颜色分量,这里将它们转化到HSV空间,获得H分量,并计算它的直方图。这里把H分量的数值量化到[0,255]。
4 实验过程
Mean Shift算法将颜色直方图作为跟踪特征,具有计算简单、速度快的特点。给定目标模板并给出t-1帧中目标的位置y0,对第t帧图像进行目标搜索的过程如下:
1)初始化第t帧中目标位置y0,并计算以及;
2)根据式(6)计算{wi}i=1,2,…,m;
3)根据式(8)计算出候选目标的更新位置ŷ1;
6)如果‖ŷ1-ŷ0‖<ε则停止,否则将ŷ1赋值给ŷ1,转至步骤2)继续。
5 实验结果与分析
为了测试本文算法在动态场景下对目标手势跟踪的性能,选择对实际视频序列进行手势跟踪测试。实验中对各帧图像未做任何预处理。计算机配置为P3 CPU 1.8GHz,内存512 Gbyte,开发环境为Visual C++。图1为实验室实时拍摄的手势跟踪结果。可以看出该算法在手势大小变化、手势发生旋转情况下均有较好的跟踪结果。
6 小结
笔者将Mean Shift方法应用到序列图像的手势跟踪中。实验结果表明,针对不同的动态场景,Mean Shift计算简单,实时性好,并且当目标姿态和光照发生变化时可获得较好的跟踪效果,能适应目标在图像中的复杂运动,具有广阔的应用前景。
摘要：将Mean Shift算法应用于序列图像中的手势跟踪,利用梯度优化方法实现快速目标定位,能够对非刚性目标实时跟踪,并且对目标的变形、旋转等运动有较好的适用性。实验结果表明,Mean Shift算法在目标姿态变化、光照变化下的跟踪效果较好。
关键词：Mean Shift算法,手势跟踪,颜色特征
参考文献
[1]WU Y,HUANG T S.Robust visual tracking by integrating multiple cues based on co-inference learning[J].International Journal of Computer Vision,2004,58(1):55-71.
[2]TRIESCH J,MALSBURG C V D.Self-organized integrationof adaptive visual cues for face tracking[C]//Proc.the Fourth International Conference on Automatic Face and Gesture Recognition.Grenoble,France:[s.n.],2000:102-107.
[3]COMANICIU D,RAMESH V,MEER P.Kernel-based object tracking[J].Pattern Analysis and Machine Intelligence,2003,25(5):564-577.
[4]COMANICIU D,RAMESH V,MEER P.Real-time tracking of nonrigid objects using mean shift[J].IEEE Computer Vision and Pattern Recognition,2000(2):142-149.
[5]CHENG Y.Mean-shift,mode seeking,and clustering[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1995,17(8):790-799.
[6]梁静,支琤,周军.基于Mean Shift的抗遮挡运动目标跟踪算法[J].电视技术,2008,32(12):82-85.
[7]KAILATH T.The divergence and Bhattacharyya distance measures in signal selection[J].IEEE Trans.Comm.Tichnology,1967(15):52-60.
[8]RICHARD O D,PETER E H,DAVID G S.Pattern Classification[M].2nd Ed.李宏东,姚天翔,译.北京:机械工业出版社,2003.
医学序列图像 篇5
关键词：离散混沌系统,连续混沌系统,图像加密,信息安全
1. 引言
自1997年J.Fridrich首次提出利用混沌加密图像[1]。利用混沌加密图像成为图像加密的一个热点,先后出现了利用Logistic映射、改进的Logistic映射、Henon映射和Chebyshev映射等离散混沌系统和lorenz、Chen、Lü等连续混沌系统,以及这些系统与M序列、传统加密方法的组合来产生用于加密图像的密钥流[2,3,4,5,6,7,8,9,10,11]。其中,Logistic映射、Henon映射,特别是Logistic映射是目前研究最广泛的混沌系统。随着研究的深入,利用单一Logistic映射加密的系统安全性受到了质疑[12,13,14],这主要是由于Logistic映射属于低维混沌系统,存在分布不均匀以及现有混沌序列发生器存在有限精度引起的退化为周期信号等原因引起的[4,15,16]。因此,寻找高维、性态复杂、分布均匀的混沌系统,从而获得大周期或非周期密钥流成为混沌加密图像的一个研究新动向。目前,大多采用多种混沌组合、超混沌、常规加密方法与混沌多级加密的方法来获得大周期或非周期密钥流。如文献17构造了一个超混沌Henon映射用于加密图=像;文献18将Lorenz、Chen、Lü三个三维连续混沌系统组合生成密钥流;文献19给出了利用Logistic、Cubic、Cat三个映射的组合来加密图像的例子;文献20则给出了将Logistic映射和Chebyshev映射组合形成密钥流的介绍。
混沌加密方法属于对称加密体制的范畴。这种加密体制的安全性取决于密钥流发生器所产生的信号与随机数的近似程度[21],密钥流越接近随机数,安全性越高,反之则容易被破译。本文介绍了一种基于离散和连续的混合混沌系统的图像加密方法,利用该方法加密的保密性能主要取决于混沌序列本身而不是传统的运算复杂度。该方法与现代密码体制的要求是一致的,即系统的保密性不依赖于对加密、解密算法和系统的保密,而仅仅依赖于密钥的保密性[22],这是本方法的最大特点。
2. 混合混沌序列的生成
混合混沌序列的生成模型如图1所示:
将两个离散混沌系统产生的混沌信号经一定的比例混合得到一个驱y'=动信y号+p(t),然后y'=再利y用+p(t)驱动连续混沌系统,从而得到系统的初始值,通过不断的迭代,得到一个新的混沌序列。该混沌序列在相图上与原连续混沌系统相似,但实际却又是完全不同的两个混沌序列,这样可以抵御基于相空间重构的攻击。
本文中,离散混沌系统I、II分别选用二维超混沌Logistic映射和Henon映射,连续混沌系统选用Lorenz系统。其系统方程如下:
1.超混沌Logistic映射:
2.Henon映射:
3.Lorenz系统:
上面各式中参数如下:图1中,
其中,L x(t)为Logistic系统,H x(t)为Henon映射,h为混合比例因子。参数h的选择以生成新的序列是否仍具有混沌性质,有无周期性为原则。本文中,经过实验验证,h取1.0,既可以满足新的序列呈现混沌特性,又能使新系统的相图和单纯的Lorenz系统相似,从而增加了破译的难度。
将式(4)代入连续混沌系统(3)得:
其中,x'=x+p(t),y'=y+p(t),z'=z+p(t)。式(5)中,
当二维Logistic映射的初值取x=00.2 y=,00.3,Henon映射的初值取x=00.65,y=00.85,Lorenz系统的初值分别取(x1,y1,z1)(=-15.0,-10.0,3.0),得到的数值模拟结果如图2所示。图2(a)为未加入驱动y'信=号y平+p(t)时,x'=方x向+的p(频t)谱图。图2(b)为加入驱y'动=信y号+p(t)时,x'=方x向+的p(频t)谱图。由图示可知,利用本方法产生的混沌序列非常接近白噪声的宽谱特性,这是本文采用Logistic映射、Henon映射和Lorenz系统组合的重要原因。限于篇幅,Logistic映射、Henon映射的频谱图没有给出,实际上它比Lorenz系统的频谱还要窄。图2(c)为加入驱动y'=信y号+p(t)时,新系统的x-z方向的相图。必须指出的是,改变混合比例因子的大小可改变该系统的相h图形状,甚至退化为周期序列,如图2(d)所示,混合比例因子h=5.0时的相图。
3. 混合混沌序列在图像加解密中的应用
由混沌序列转换为加密因子的原则是让加密因子序列尽量保持原混沌序列的伪随机性,同时适于加密算法的操作。目前,该变换大体上可以分为实数值序列、位序列和二值序列三种[17]。本文采用的方法是将混沌实数序列按256个等级线性等分,然后映射得到8位二进制序列,将该序列作为加密因子序列。对图像进行加密时,将该加密因子序列,依次对图像的每一个像素值进行异或运算(模2加)便可得到加密图像。这种算法简单、计算速度快,且易于实现。解密过程与之相逆。图3为利用混合混沌序列中的x方向的数值进行加密、解密的实验结果。
4. 本方法的特点
加密算法的安全性主要依赖密钥的安全性,而密钥的安全性则与密钥空间的大小,密文对密钥的敏感性以及算法的复杂性都有很大的关系[23]。下面我们通过理论分析和实验结果来评估该系统加密的安全性。
4.1 密钥空间
混合混沌系统的参数a1,a2,a3,a4,a5,a6,h,σ,ρ,b和各系统的初始值都可以包含在密钥空间中。若运算精度按照10-16考虑,那么Lorenz系统的初始条件组成的空间大小至少为1016×1016×1016=1048。Logistic二维离散超混沌映射和Henon映射的初始值又可以组成至少是1016×1016×1016=1048的空间。另外还有a1,a2,a3,a4,a5,a6,h,σ,ρ,b等参数都可以作为密钥。因此,本系统具有非常巨大的密钥空间,要进行穷尽攻击是很难的。
4.2 加密图像对初始条件的敏感性
一个好的加密算法应该实现密文对密钥的敏感性,即密钥的微小差别可以对密文造成很大的变化。密文对密钥越敏感,密文中包含的密钥信息量就越少,使得基于唯密文的攻击越困难,因为从密文的变化中很难得出密钥的变化情况。图3(c)是密钥一致时的解密结果,图3(d)是密钥值相差10-14时的解密结果,可见,加密图像对密钥是非常敏感的。
4.3 算法抵抗基于统计分析攻击的能力
直方图是反映像素分布的一个重要指标,在此,将原始图像的直方图和加密后图像的直方图进行比较。从图4(a)和图4(b)可知,加密图像的直方图具有更好的均匀特性。再来看加密图像相邻像素的相关性分析。我们通过计算原图像和加密图像水平、垂直和对角相邻像素的相关系数来分析。图4(c)给出了加密前后图像的水平方向相邻像素的相关系数。可以看出,加密后相关性大大下降。实验结果表明,该方法对图像具有较好的扩散和混淆能力,能抵抗统计分析的攻击。
5. 结论
医学序列图像 篇6
当今的网络由于其协议、系统方面的缺陷, 一些有特定意义的图像 (如医疗图像、法庭证据等) 易受到非法者的篡改、盗取, 图像加密算法能有效保护图像拥有者的利益。
混沌系统广泛应用在图像加密系统中[1,2,3], 但文献[4]表明单一的混沌映射图像加密算法安全性较差。
Adleman通过生物实验成功求解7个顶点的有向汉密尔顿路径问题[5], DNA计算诞生。Catherine等人[6]提出一种以DNA为载体来隐藏信息的方法。Gehani等人[7]提出一种基于DNA序列特性的一次一密加密算法。此类方法需要复杂的生物实验来实现, DNA密码学多停留在理论层面。Kang Ning等人[8]提出一种伪DNA计算的加密方法, 不需要复杂的生物实验, 还具有很好的加密效果, 自此提出了诸多基于DNA特性的加密算法[9,10,11]。
本文算法, Chebyshev混沌对原图进行基于位平面的像素置乱[1], 利用天然DNA序列作为密钥, 根据DNA的互补配对原则, 自定义DNA矩阵的加、减运算, 再结合PWLCM混沌映射进行像素替代。
1. 相关理论
1.1DNA编码规则
天然DNA序列由4种碱基A、T、C、G组成。根据碱基互补的原则, 定义8种DNA编码规则, 如表1。
8bit的灰度图像, 一个像素点可由4个碱基来表示灰度值。
1.2DNA矩阵运算
表2、表3示例采用表1中第一种DNA编码规则定义一种DNA序列的加、减运算:
1.3混沌映射
Chebyshev混沌映射:
α>2, xn∈[-1, 1], 初始值变化为1*10-16时可产生两个截然不同的混沌序列。
PWLCM混沌映射:
xn∈[0, 1) , p∈ (0, 0.5) , 该混沌具有均匀的不变分布以及更高的遍历性、混合性和确定性。
1.4混沌映射初始值获取算法
首先将原图每个像素点转换为8位的2进制序列, 再划分为4个坐标对, 定义集合{A}、{B}一一对应存储坐标对的值。例如图像像素点为54、92……, 集合{A}、{B}对应值如图1。再计算得到新的集合{a}、{b}, 其中a0=b0=0.5。
最后计算得到混沌映射的初始值x0和y0, 其中m、n是图像的行列维数。
算法中通过PWLCM混沌映射的序列对得到的DNA序列进行碱基互补, 混沌映射的初始值是由天然DNA序列根据上述原理获得, 使得算法的像素替代敏感于天然DNA密钥。
2. 算法流程
加密算法流程如图2所示。
1) 读取灰度图像, m、n是图像I的行列, 根据1.4节的方法得到Chebyshev映射的初始值x0和y0。
2) Chebyshev映射采用初始值x0迭代m次, 得到序列{X}, 然后对序列{X}进行升序排列得到{X'}, 将图像I扩展后的位平面的行根据序列{X}各元素在{X'}中的位置进行排序。
Chebyshev映射采用初始值y0迭代n*8次, 得到序列{Y}, 然后对序列{Y}进行升序排列得到{Y'}, 将图像I扩展后的位平面的列根据序列{Y}各元素在{Y'}中的位置进行排序。
最后得到置乱后的图像A (m, n) 。
3) 置乱图像A (m, n) 扩展为8个位平面, 两个相邻位平面为一组, 每一组按照1.1节定义的原则分别采用1~4种规则进行DNA编码, 得到4个DNA矩阵A1 (m, n) , A2 (m, n) , A3 (m, n) , A4 (m, n) 。
4) 天然DNA序列, 生成4个天然DNA矩阵NM1 (m, n) , NM2 (m, n) , NM3 (m, n) , NM4 (m, n) 。
5) 根据1.2节定义的DNA序列加法, 矩阵A1, A2, A3, A4和NM1, NM2, NM3, NM4各自对应相加得到矩阵C1, C2, C3, C4。
6) 矩阵NM1与NM2拼接得到矩阵N1 (m, n*2) , 矩阵NM3与NM4拼接得到矩阵N2 (m, n*2) , 矩阵M1和M2按照1.1节的第一种编码规则进行解码得到二进制矩阵M1 (m, n*4) 和M2 (m, n*4) , 然后M1和M2按照1.4节的原理得到初始值a0, b0, c0, d0。
7) PWLCM映射, 采用初始值a0, b0, c0, d0, 分别迭代m*n次, 得到4个序列{P1}, {P2}, {P3}, {P4}, 然后将4个序列转为二值序列, 当pij=1, i∈[1, 2, 3, 4], j∈[1…mn]时, DNA矩阵Ci对应位置的碱基取补, 否则不变。得到矩阵C1’, C1’, C1’, C1’。
8) 最后将C1’, C1’, C1’, C1’按照1.1节定义的原则分别采用5~8种规则进行DNA解码, 再将解码得到的二进制矩阵依次按照8个位平面组合, 最终得到加密的灰度图像EI (m, n) 。
解密算法是加密算法的逆操作, 步骤5的加法采用减法。
3. 实验结果及分析
在matlab7环境下, 对3幅128*128像素的灰度图像“Lena”、“Babara”、“Peppers”进行仿真实验。密钥key={AL031432, 1, 4, 0.35}, 天然DNA序列Gene ID是AL031432, 从第1位碱基截取128*128*4的DNA序列, 4是Chebyshev混沌映射参数α, 0.35是PWLCM混沌映射参数p。图3是加密图像。
3.1置乱程度分析
1) 不动点比
原图I (m, n) , 加密图I' (m, n) , 公式:
2) 灰度平均值
原图I (m, n) 和加密图I' (m, n) , 灰度平均值计算公式如下:
灰度平均值产生均匀变化时, 安全性高, 置乱效果好。
3) 信息熵
信息熵是测试图像灰度值分布的一个参数, 公式如下:
256灰度级的灰度图像, 如果每一个灰度值出现概率都相同, 信息熵H (m) =8, 即该图完全随机。为了抵御基于熵的攻击, 加密图像信息熵越接近8越好。
进行仿真实验, 表4给出了置乱程度定量分析的结果, 结果表明该算法置乱效果非常好。
3.2相关性分析
一个有效的加密算法能降低相邻像素点的相关性。随机选取4000对相邻像素点 (横向、纵向、对角线) 进行相关性分析。分析函数:
表5是图像加密前后的相关性信息。本文的算法能有效抵御基于相关性分析的攻击。
3.3差分攻击分析
NPCR (像素变化率) 和UACI (归一化平均变化强度) 这两个标准来衡量加密算法是否能有效抵御差分攻击, 公式:
原图改变一个像素点Lena:T (1, 1) =255、Babara:T (3, 3) =10、Peppers:T (2, 2) =250, 加密后进行分析, 表7, 由结果可知, 本文提出的算法能有效抵御差分攻击。
3.4密钥空间及敏感性分析
一个好的加密算法必须有极大的密钥空间才能有效抵御穷举攻击, 随着首个完整基因组测序工程的完成, DNA序列数据呈现指数级的增长, 且对同一个DNA序列, 截取起始位置不同, 获取的序列完全不一样。
加密算法对密钥敏感性越高该算法的安全性越高。仿真测试分别改变密钥DNA矩阵的第一位碱基、第一百位碱基、第二百位碱基为其互补的碱基, 然后对各图进行错误密钥解密, 图4展示了错误密钥解密后的图, 与原图完全不同, 表示本算法对密钥有极高的敏感性。
4. 总结
本文提出的基于DNA序列和混沌映射的图像加密算法, 通过原图获取Chebyshev混沌初始值, 加密算法敏感于原图, 基于位平面置乱图像像素点, 实验结果表明置乱程度高;天然DNA序列作为密钥, 密钥空间大, 能模拟一次一密的加密效果, 且密钥的存储、传输更为便捷;置乱后图按位平面划分生成DNA矩阵与天然DNA矩阵运算, 且各平面的编码规则各不相同, 算法加密性能好;最后通过PWLCM混沌序列, 对得到的DNA序列进行碱基互补配对, 配对操作敏感于天然DNA序列, 算法安全性高。实验结果表明, 本文提出的加密算法, 置乱程度高, 安全性高, 密钥空间大, 能有效的抵御穷举攻击、统计攻击、差分攻击。
参考文献
[1]Fu C, Lin B, Miao Y, et al.A novel chaos-based bitlevel permutation scheme for digital image encryption[J].Optics Communications, 2011, 284 (23) :5415-5423.
[2]Awad A, Awad D.Efficient image chaotic encryption algorithm with no propagation error[J].ETRI journal, 2010, 32 (5) :774-783.
[3]Zhu C.A novel image encryption scheme based on improved hyperchaotic sequences[J].Optics communications, 2012, 285 (1) :29-37.
[4]Solak E, Rhouma R, Belghith S.Cryptanalysis of a multichaotic systems based image cryptosystem[J].Optics Communicatio ns, 2010, 283 (2) :232-236.
[5]Adleman L M.Molecular computation of solutions to combinatorial problems[J].SCIENCE, 1994:1021-1021.
[6]Clelland C T, Risca V, Bancroft C.Hiding messages in DNA microdots[J].Nature, 1999, 399 (6736) :533-534.
[7]Gehani A, La Bean T, Reif J.DNA-based cryptography[M].In Aspects of Molecular Computing.Springer Berlin Heidelberg, 2004:167-188.
[8]Kang Ning.A pseudo DNA cryptography method[J].DBLP:journals/corr/abs-0903-2693, 2009.
[9]Zhang Q, Guo L, Xue X, et al.An image encryption algorithm based on DNA sequence addition operation[C]//Bio-Inspired Computing, 2009.BIC-TA'09.Fourth International Conference on.Ieee, 2009:1-5.
[10]Zhou Shihua.The Research of Image Encryption Based on DNA Sequences and Self-assembly[D].Dalian University of Technology, 2013.
医学序列图像 篇7
关键词：图像分割,二维直方图加权模糊聚类算法
1、引言
图像分割在乳腺MR图像的定量分析和定性分析中扮演着十分重要的角色, 直接影响后续处理与分析工作[1,2]。通常动态增强图像的分析依赖于对感兴趣区域的分析, 感兴趣区域的选择有多种方法[3,4]。针对乳腺MRI图像病变分割的难度, 多种分割方法得到广泛的应用[5,6,7,8]。本文提出一种基于二维直方图的加权模糊C-均值聚类算法分割肿瘤, 该方法实际上对增强减影序列图进行分割, 而不是直接在时间序列图上进行, 主要是因为减影图能够保证增强区域与正常组织之间的对比度, 该方法能够快速地提取肿瘤的轮廓。
2、DWFCM算法原理
传统FCM算法基于灰度信息的分类算法, 仅利用了图像的灰度信息, 然而灰度信息表达有限, 尤其对于模糊MRI医学图像, 灰度级之间的区分度较难, 这使得传统FCM方法并不能达到很好的分割效果, 因此本文提出2DWFCM肿瘤分割方法, 该方法不仅考虑了图像像素点的灰度信息和其邻域的均值信息, 同时还考虑到不同样本点对于分类的贡献不同来赋予适当的权重, 从而达到较好的肿瘤分割效果。
2.1 二维直方图
由于每一图像的像素与其邻域像素的相关性相当大, 当利用像素空间相关的信息, 通过对原图进行邻域平均得到另一幅平滑的图像, 由原图像和平滑图像构造出一个二维的直方图, 物体和背景的分布在二维直方图中就会比在一维直方图更容易区分。因此, 利用图像的空间信息进行分割一定会产生更好的效果。
设图像的尺寸大小为M×N, 图像灰度级数为L, 显然图像中坐标 (m, n) 像素点的灰度值I (m, n) ∈[0, L-1]。定义坐标 (m, n) 处像元点的邻域灰度均值g (m, n) 如下:
这里S为正方形邻域窗口的宽度, 一般取奇数, 本文中取3。其中0<m<M, 0<n<N且g (m, ) nD∈[0, L-1]。对于任意一帧图像T=T (m, n) , 由I (m, n) 和g (m, n) 可以构成一个二元组, 可以用二维向量矩阵表示:
对于一幅MXN的图像, 采用以上二维向量表示时, 定义并计算它的二维直方图。该二维直方图定义在一个L×L大小的正方形区域内, 二维坐标分别表示图像像元的灰度值和邻域灰度均值。直方图中任意一点的值定义为pij, 表示向量 (i, j) 发生的概率, 其中 (i, j) 表示 (I (m, n) , g (m, n) ) .用hij表示向量 (I (m, n) , g (m, n) ) 发生的频数, 那么向量 (I (m, n) , g (m, n) ) 发生的概率表示如下:
其中, 并且显然有如下关系式:
2.2 FCM聚类算法
FCM是一种无监督的机器学习方法, 它的目标是将具有N个点的数据集X模糊划分成C类。若将图像像素组成n个样本集合X{xi, , i=1, 2…, n}, 那么图像的分割问题就转化为将这n个样本分成C个聚类的问题。Vk{k=1, 2, …, n}为每个聚类的中心, uik是第i个样本对第k类的隶属度函数。J.C.Bezdek提出的FCM聚类算法是一种模糊目标函数法, 其目标函数J (U, V, X) 定义为:
其中, m∈[1, +∞) 是加权指数, 用来控制分割结果的模糊性, ||·||代表欧氏距离, U=[uik]则为隶属度矩阵, 满足如下关系式:
J (U, V, X) 表示各类别样本到聚类中心的加权误差平方和, 其值越小, 聚类越紧凑, 这实际上是一个数学规划问题, 利用拉格朗日乘数法对 (2-5) 式求条件极值, 得出以下两个公式:
其中, k=1, 2, …, i=1, 2, …, n, 。
FCM算法已被证明具有良好的收敛性, 收敛后的就是实现二值化:假设背景的平均灰度值小于背景的, 那么当||v1||<||v2||时, v1和v2对应背景和物体的聚类中心;反之, 则v1和v2对应物体和背景的聚类中心。对于第一种情况的二值化过程如下:%
其中, α为二值化的阈值, 可选为[0.5, 1]间的某一个值, uik≥α表示第i个像素点隶属于第k类的程度大于等于α。
2.3 2DWFCM算法
基于二维直方图加权的模糊C均值算法首先考虑到不同像素点对分类的影响, 给予每个像素点一个权重ωi{i=1, 2, …, n}, 根据这一原理, 将目标函数 (2-5) 修正如下:
利用拉格朗日乘数法对 (2-10) 式求条件极值, 得出以下迭代公式:
其中k=1, 2, …, i=1, 2, …, n, ωi是每个像素点所对应的加权系数, 满足概率约束条件:
由此得出, 隶属度的迭代公式与经典的FCM中的隶属度迭代公式是一致的, 不同的是聚类中心vk的迭代公式, 因此, 权重系数ωi主要作用于聚类中心的调整。我们定义待分类的样本为T=[I (m, n) , g (m, n) ], 加权系数为ωi=pij, 也就是将对应点的二维直方图的值作为对应样本的权重。
3、肿瘤分割
本文利用2DWFCM算法从乳腺DCE-MRI图像中提取肿瘤, 肿瘤分割过程主要包括六个连续的步骤:获取动态增强图像的剪影图像;用户选择感兴趣区域;利用2DWFCM算法提取肿瘤;肿瘤二值化;连通标记及目标选择;目标空洞的填补;图1即为该分割的程序框图。这里我们用一个实例来详细说明分割的步骤.
如图2实例图可知每个时间序列包括多张序列图, 其中上方的三幅图为第一张出现病灶、出现病灶最明显和最后一张出现病灶的图片, 矩形框标出的为病灶的位置。整个分割过程用户的干涉只涉足到感兴趣区域的选择, 用户通过选择对话框来选择病灶最明显的切片, 在此切片上拖动鼠标截取一个包含病灶的矩形子图。对于一个完整的乳腺DCE-MRI序列图实际上用户只需要操作一次, 其他的序列图复制该位置即可, 其中下方为增强前及注射增强剂后六个不同时间点采集的含有病灶的序列图。
乳腺病变组织较正常组织而言, DCE-MRI强化很明显, 为了突出病灶, 用增强后的序列图减去相应增强前的序列图, 得到的结果序列图作为后续处理的数据源。
利用2DWFCM算法对增强的时间减影序列进行分类, 将其分为两类 (C=2) :病灶与非病灶, 分类的实质就是最小化2-10式的目标函数。首先随机地初始化聚类中心矩阵V, 隶属度矩阵U初始化为0。然后利用式子2-11和2-12迭代更新矩阵U和V直至算法收敛。这里设置收敛参数为一个很小的数ε (ε=10-6) , 加权指数m设为m=2。2DWFCM的分割结果图如3所示。该图分别显示了肿瘤病灶提取的各个步骤的结果:利用2DWFCM分割后选择一个合适的阈值进行二值化, 根据经验这里阈值取0.2较好, 如图3 (a) 所示。由图3 (a) 可以发现肿瘤病灶的周围存在一些假阳性点, 为了剔除这些假阳性点, 本文对3 (a) 的结果图进行连通标记操作及目标的选择, 实现过程如下: (1) 扫描二值图像的所有像素点, 为所有非零点初始化标签。 (2) 对 (1) 中非零点划分等价类。 (3) 根据 (2) 中得到的等价类重新给每个点贴标签。 (4) 对 (3) 中标记的连通区域进行目标的选择, 在这里选择包含病灶中心的连通区域作为肿瘤病灶组织, 如图3 (b) 所示, 其他所有连通区域像素点的结果赋值为零。目标选择完成后将执行漏洞的填补操作, 这主要是因为在肿瘤病灶区域可能存在一些坏死的组织, 这些坏死的组织在增强过程中未得到增强, 因此最初在用2DWFCM算法进行分割的时候这些区域就被分成非病灶类。为此采用形态学方法来填补这些漏洞, 如图3 (c所示。
4、实验分析与评价
实验数据源由福建省肿瘤医院提供, 这里用到的图片为病人增强前的序列图与增强后的时间序列图。病人的序列图为三维T2加权梯度回波序列图, 相关的参数为:重复时间TR=5.9 ms, 回波时间TE=2.8 ms.在获取增强前扫描序列后, 以2m L/s的速度给病人注射02mmol/kg的Gd-DTPA, 随后以同样的速度注射20m L盐水冲洗。注射增强剂后以60s的时间间隔连续扫描六个序列。每个序列包括62张乳腺切片, 每张切片图的大小为512×512。切片的厚度范围为2.0~3.0mm, 因乳腺的大小而定。
图4为实验中部分病人的肿瘤分割结果图, 其中第一行为8个病人的感兴趣区域, 其中四个为良性肿瘤感兴趣区域, 如Case1~Case4;五个为恶性肿瘤感兴趣区域, 如Case5-Case8。第二行为采用本文提出的方法2DWFCM对这8个病人的感兴趣区域分割的相应结果图, 由图可知本文的方法无论是对良性肿瘤还是恶性肿瘤的分割都能获得较好的效果。第三行为采用FCM方法对这8个病人的感兴趣区域分割的相应结果图, 这些分割结果与2DWFCM分割的结果相比, 有的错将周边的假阳性划为肿瘤的范畴, 有的则产生了过分割, 尤其是对浸润性的恶性乳腺肿瘤, 如图Case7分割图所示。这主要是因为浸润性的恶性肿瘤内部的内部结构不均一, 常常跟周围正常组织浸润在一起, 甚至有很多坏死的组织, 因此增强不明显, 导致FCM方法分割效果不是很好。由FCM与2DWFCM的分割结果对比表明, 本文的方法对动态增强的乳腺MRI图像肿瘤的分割具有较好的分割效果。
5、结论
本文针对对比动态增强乳腺MRI图像的特点, 提出基于二维直方图的加权模糊聚类的分类方法, 将含有肿瘤的感兴趣图划分为两类:肿瘤与周围正常组织。并且将实验结果与FCM方法的实验结果进行了对比, 实验结果表明了该方法的有效性。本方法的缺陷是未能够达到完全自动化的过程, 今后的改进方向是寻找如何自动的获取感兴趣区域, 而无需人工干预。
参考文献
[1]田捷, 包尚联, 周明全.医学影像处理与分析.电子工业出版社.2003.
[2]Chen W, Giger ML, Lan L, Bick U.Computerized interpreta-tion of breast MRI:investigation of enhancement-variance dynam-ics.Med Phys 2004;31:1076-82.
[3]Liney GP, Gibbs P, Hayes C, Leach MO, Turnbull LW.Dy-namic contrast-enhanced MRI in the differentiation of breast tumors:user-de?ned versus semi-automated region-of-interest anal-ysis.JMagn Reson 1999;10:945-9.
[4]ChenW, GigerM, Bick U, Newstead G.Automatic identi?cationand classi?cation of characteristic kinetic curves of breast lesions onDCE-MRI.Med Phys 2006;33 (8) :2878-87.
[5]Chen W, Giger ML, Bick U.A fuzzy c-means (FCM) -basedapproach for computerized segmentation of breast lesions in dy-namic contrast-enhanced MR images.Acad Radiol 2006;13:63-72.
[6]Coto E, Grimm S, Bruckner S, Groeller E, Kanitsarand A, Ro-driguez O.MammoExplorer:an advanced CAD application forbreast DCE-MRI.In:Proceedings of Vision, Modelling, and Visu-alization.2005.
[7]Lucht R, Delorme S, Brix G.Neural network-based segmenta-tion of dynamic MR mammographic images.Magn Reson Imaging2002;20:147-54.