本发明属于医学图像处理，具体涉及一种基于em算法与贪心策略的红细胞序列样本分割方法。
背景技术：
科学技术的发展使得对红细胞的研究越发向微观方向迈进，且研究方法也越发多样。计算机的发展更是极大的推动了研究的进步。在此基础上，科研者们进行了诸多的研究与探索：何晓川，崔建强等人对红细胞的流动过程进行了基础性的研究[1][2]，他们建立了红细胞的模型，分析了红细胞与相邻流体的相互作用，最后利用在流体力学领域广泛应用的软件fluent对红细胞的流动过程进行模拟，并着力于分析形成的流场。在此基础上，杨婷[3]探究出了一种研究红细胞形变的新方法：她将红细胞模拟为含有牛顿流体的双凹面碟形弹性薄膜胶囊，然后采取多块策略改善细胞模型附近的网格，提高了三维计算的效率及结果的准确度。陈春晖[4]总结了前人的研究成果，在此基础上加入水平集法追踪仿真血液中的红细胞。结果表明这种方法可模拟红细胞的运动，且准确性也得到了极大的提升。随后，赵立国将任意拉格朗日方法(ale)和levelset方法相结合[5]：引入ale方法处理仿真中的边界大尺度畸变问题，在保证计算精度的同时增强了灵活性；levelset方法则是将低维闭合曲线中的问题转化到高维空间中，从而追踪模拟活动边界。结果表明，该方法可应用于复杂拓扑结构变化。([1]何晓川.微循环中红细胞流动的基础性研究[d].济南：山东大学，2009.[2]崔建强.微血管红细胞流动特性的理论与仿真[d].济南：山东大学，2010.[3]杨婷.红细胞在微小血管中的流动理论分析[d].济南：山东大学，2011.[4]陈春晖.毛细血管中红细胞流变特性的基础性研究[d].济南：山东大学，2012.[5]赵立国.红细胞在毛细血管中的流动形态分析与仿真[d].济南：山东大学，2014.)
目前对红细胞流动的仿真模拟方法还是侧重于在流体力学的基础上通过不同的方法对血浆运动进行模拟。
尽管技术的发展使得红细胞相关的研究与仿真水平相较于从前已经得到了很大程度上的提升。但仍存在诸多不足：对细胞仿真研究现状的分析显示，以往的仿真方法多集中于对各类细胞静态模型的建立。但红细胞的形态会因为外部条件的变化而产生不同的形变，静态模型方法难以展现这类变化的过程，所以需要实现动态过程的仿真。对红细胞流动形态仿真研究现状的分析显示，目前对红细胞流动的仿真重点是在流体力学的基础上通过不同方法对血浆运动进行模拟，而我们希望能既展示出大量红细胞的整体流动过程，又展示出在不同的外界条件下单个红细胞形态的变化过程。因此，引入动态纹理技术，以实现对红细胞流动过程的更具真实性的仿真。
技术实现要素：
本发明的主要目的在于提供一种基于em算法与贪心策略的红细胞序列样本分割方法提高基础医学仿真的真实性与直观性，在此基础上减少仿真过程中的手动工作量。
为实现上述目的，本发明采用以下技术方案：一种基于em算法与贪心策略的红细胞序列样本分割方法，其包括以下步骤：
步骤s1：将红细胞序列样本中所有可直接观测到的状态集合表示为y，不可直接观测到的状态集合表示为x，用θ作为描述所有状态集合联合概率的模型参数，并以p(x,y|θ)对两状态集合的联合概率进行表示；以y表示观测所得数据，则显然y是已知量；模型参数利用最大似然估计法实现估计，如式(1)所示：
步骤s2：对式(1)展开等价变化得到基于初始分布
1.e步：更新不可直接观测到状态的分布，即：
2.m步：更新模型参数
其中以t＝1,2，…表示迭代次数；
步骤s4：取一定长度tmin的红细胞流动样本帧序列对步骤s3中形成的em算法形成线性系统进行训练；
步骤s5：判断可靠性指标r是否小于给定值，若是添加后续帧，否则执行步骤s6
步骤s6：将此线性系统在其余未标记的序列上进行测试，看是否有红细胞流动样本帧序列与之相对应，即红细胞流动样本帧序列是否符合同一个线性系统，视频中包含同一纹理基元的多个复制；若没有，则需利用添加该帧序列后的新序列训练一个新的线性系统；
步骤s7：对以上步骤s4-s6不断地重复，直到将样本序列分割完毕。
在本发明一实施例中，步骤s2包括以下具体步骤：对式(1)展开等价变化如式(2)所示：
其中ep(x|y,θ){.}用于表示随机变量x在分布为p(x|y,θ)的期望，但实际上，x的真实分布密度x～pθ＝p(x|y,θ)属于未知量，所以基于初始分布
其中d()为概率分布函数、h()为概率密度函数。
与现有技术相比，本发明具有以下优点：
1.在红细胞的形态仿真过程中，相似的动态纹理基元可能会多次重复出现，对其的分割与重新拼接可能造成对原始运动的精确重复，这就会使合成结果缺乏真实感与自然感。本发明最终合成结果应避免此问题的出现。
2.若分割出的纹理基元并非位于序列中相邻位置，那么直接跳转会使合成结果中的转换过渡显得生硬，缺乏真实感。但若对其进行平滑滤波，又会破坏流动过程的一致性。本发明最终合成结果应实现真实感与一致性的统一。
3.对于拓展过程中不可避免的噪声问题，最终合成结果应实现噪声影响的最小化，即能在一定程度上减少噪声累计，从而避免失真现象。
附图说明
图1为本发明一实施例样本分割过程流程图。
具体实施方式
下面结合附图和具体实施例对本发明做进一步解释说明。
对于正常状态下的红细胞流动仿真序列，由于其本身具有规律性，因此只需要较短的样本就能够合成无限长且无视觉重复的动态序列。而对于外界条件改变的情况，由于规律性较弱，因此需要较长的样本序列，且合成的序列是有一定长度的。
本发明采用融合了基于em算法与贪心策略的方法实现对红细胞序列的样本分割。
em算法，即最大期望算法，是一种迭代算法。采用em算法对不可直接观测到的状态模型进行估算。em算法的核心思想是在参数未知的情况下假定一个参数初值，然后根据这个初值进行不断迭代，对模型参数进行不断更新直至其收敛。
em算法的每一个迭代周期包含以下两个步骤：
1.e(expectation)步骤：主要是用于计算期望。假设当前模型参数无误，采取统计推断算法，结合观测数据计算隐节点的分布。
2.m(maximization)步：主要用于实现最大化。假定r步得到的统计分布是正确的，可以直接利用最大似然准则更新模型的参数。
em算法的具体操作如下：若将样本序列中所有可直接观测到的状态集合表示为y，不可直接观测到的状态集合表示为x，用θ作为描述所有状态集合联合概率的模型参数，并以p(x,y|θ)对两状态集合的联合概率进行表示。以y表示观测所得数据，则显然y是已知量。模型参数可以利用最大似然估计法实现估计，如式(1)所示。
对式(1)展开等价变化如式(2)所示。
以ep(x|y,θ){.}用于表示随机变量x在分布为p(x|y,θ)的期望，但实际上，x的真实分布密度x～pθ＝p(x|y,θ)属于未知量，所以基于初始分布
用f(p,θ)表示e,m两步迭代过程，即有
1.e步：更新不可直接观测到状态的分布，即
2.m步：更新模型参数
这里以t＝1,2，…表示迭代次数。
在实际操作当中，若使每一步都达到最优化，不但会造成时间的消耗，而且会使得实现难度增加，在可操作性方面并不乐观。因此，在操作之中可不需使em迭代过程的每一步均为最优参数值，只要确保每次迭代均可使能量函数的数值不断增大，算法就是收敛的。
基于这一点，将em算法与贪心策略相融合，实现所需效果。
贪心策略的设计特点是算法逐步进行，基于当下情形，参照某一特定优化测度做出选择，而并不对各式各样可能的整体情形进行考量。其结构为自顶向下，采取迭代的方法进行选择，每次选择都能够使得问题规模缩小，从而达到简化问题的目的。每一步贪心选择均得到问题的一个最优解。使用这一算法主要是为了避免时间上的消耗。该算法不需要寻找最优解，因此也不需要穷尽所有可能。
融合贪心策略后的具体操作如下：
1.取一定长度的红细胞流动样本帧序列对线性系统进行训练。一般来说，这一长度即为前文中所提及的tmin。
2.添加后续帧直至图像取得一定的可靠性。这一项通常以可靠性指标p是否小于给定值来进行判断。
3.将此系统在其余未标记的序列上进行测试，看是否有序列可以与之相对应，即序列符合同一个线性系统，视频中包含同一纹理基元的多个复制。若没有，则需利用添加帧后的新序列训练一个新的线性系统。
4.对以上步骤进行不断地重复，直到将样本序列分割完毕。
此过程的流程大致如图1所示。
本发明具体包括以下步骤：
步骤s1：将红细胞序列样本中所有可直接观测到的状态集合表示为y，不可直接观测到的状态集合表示为x，用θ作为描述所有状态集合联合概率的模型参数，并以p(x,y|θ)对两状态集合的联合概率进行表示；以y表示观测所得数据，则显然y是已知量；模型参数利用最大似然估计法实现估计，如式(1)所示：
步骤s2：对式(1)展开等价变化得到基于初始分布
步骤s3：用f(p，θ)表示em算法的e，m两步迭代过程，即有
1.e步：更新不可直接观测到状态的分布，即：
2.m步：更新模型参数
其中以t＝1,2，…表示迭代次数；
步骤s4：取一定长度tmin的红细胞流动样本帧序列对步骤s3中形成的em算法形成线性系统进行训练；
步骤s5：判断可靠性指标r是否小于给定值，若是添加后续帧，否则执行步骤s6
步骤s6：将此线性系统在其余未标记的序列上进行测试，看是否有红细胞流动样本帧序列与之相对应，即红细胞流动样本帧序列是否符合同一个线性系统，视频中包含同一纹理基元的多个复制；若没有，则需利用添加该帧序列后的新序列训练一个新的线性系统；
步骤s7：对以上步骤s4-s6不断地重复，直到将样本序列分割完毕。
在本发明一实施例中，步骤s2包括以下具体步骤：对式(1)展开等价变化如式(2)所示：
其中ep(x|y,θ){.}用于表示随机变量x在分布为p(x|y,θ)的期望，但实际上，x的真实分布密度x～pθ＝p(x|y,θ)属于未知量，所以基于初始分布
其中d()为概率分布函数、h()为概率密度函数。
以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。