An Extraction Algorithm for Motion Parameters of A Laboratory Mouse by Model Matching and Particle Filtering
-
摘要: 实验小鼠是一种变形体对象，现有方法难以从连续视频图像中同时提取出运动轨迹和体态细节.本文采用模板匹配和粒子滤波的目标跟踪方法求解这一问题.提出了一种几何体部件模型，在引入小鼠移动速率的基础上给出了其运动状态方程，以二值化前景像素与几何部件模型间的差异度方程为观测模型，以状态方程及相互独立的多维随机变量为运动模型，从而建立起基本粒子滤波算法.与逐帧差分识别方法的对比实验研究表明，所提出的模型与实验小鼠形体相似，能够达到视频在线提取的计算效率.新方法在强噪声干扰条件下解决了运动轨迹和体态同时精确估计，并有效避免了首尾识别混淆及虚影干扰等困境，从而为后续生物学行为分析提供依据.Abstract: Laboratory mouse is a kind of deformable object. Existing methods can hardly extract motion trajectories and posture details simultaneously from those continuous recorded videos. An object tracking method based on model matching and particle filtering is adopted to solve this problem. A geometry based part model and its motion state function involving moving velocity are proposed. A model-observation difference function is established as the observation model by comparing the foreground pixels in the binary image and the geometry part model. A basic particle filter is built with this observation function and the motion state function with multi-stochastic variables which follow an independent distribution. Comparison is made between the proposed method and the classical frame-differencing method, which proves that the novel part model is analogous with a physical mouse in shape and supports real-time extracting rate and high computing efficiency. The novel method is able to estimate precisely both motion trajectories and posture states, and avoid effectively the faults of head-tail confusion and reflection disturbance. Therefore the novel method provides a trust worthy means for later behavioral analysis for biologists.
-
Key words:
- Object tracking /
- particle filter /
- part model /
- laboratory mouse /
- posture
-
小鼠是生物学和医学研究中常用的一种实验动物.对实验小鼠的运动轨迹和体态细节跟踪能够实现其行为学分析, 对神经生物学、药理学和心理学等学科具有重要的意义[1].在行为记录和定量分析过程中, 视频监视和数字图像处理是一种应用极为广泛的方法.这是由于其具有布置简单、成本低廉、可用信息量大及非介入等优点.经摄像采集、信号转换、输出记录等一系列过程可获得连续的数字图像序列.在行为学分析时, 首先从这一数字图像序列中提取出小鼠的运动状态参数, 如活动轨迹、移动速度及体态等, 从而统计出运动频率、行为节律等生物学规律, 而自动化的运动状态参数提取将大幅提升这一过程的效率和准确性[2-3].
目前已有诸多小鼠运动监控系统[4-5]以及成熟的商业软件, 如EthoVision [6]等.这些系统一般基于小鼠与实验背景的颜色或纹理差异, 采用背景学习和背景减除的方法提取小鼠的二值化图像, 计算其包围盒和重心, 从而针对每一张数字图像给出目标小鼠的位置及粗略的方位角[5], 再经逐帧位置差分获得其移动速度, 通过连接小鼠在每一帧图像中的位置, 绘制出其移动轨迹.在细致的参数标定及滤波平滑处理条件下, 这些方法能够完成多种光照条件乃至黑暗条件下的运动识别任务[7-9], 然而这种基于逐帧识别的方法易受到前景提取误差和小鼠体态变化的干扰, 产生跳跃的运动轨迹和较大移动速度估计偏差, 从而难以满足那些对动作细节要求较高或小鼠活动频率要求精确统计的情形.此外, 这些方法常常忽略小鼠体态的运动细节, 例如蹲窝、蜷缩等.例如Pistori等[10]虽然实现了多只小鼠存在互相接触条件下的相对位置跟踪, 引入粒子滤波器实现了接触与非接触状态的区分以及各对象的宏观位置识别, 但难以解决小鼠体态信息的提取难题, 进而无法为深入分析小鼠行为提供更多依据.
张敏等[11]以轮廓曲线的曲率及其频谱作为特征矢量, 针对一定数量图像的聚类和学习, 实现了小鼠的修饰、伸长、曲身等体态识别; Fabrice等[12]通过几个基本几何元素及其相对运动关系建立了小鼠的理想模型, 通过物理学运动规律对该理想模型与二值化图像中的前景像素点进行匹配, 从而在非标记条件下获得了其头部、颈部、尾部等运动细节. Qing等[13]通过识别尾部、头部等特征点较为精确地提取了头部转角及体长等信息. Branson等[14]基于小鼠的样条轮廓模型及Gauss随机运动模型采用粒子滤波器实现了几个小鼠的轮廓最佳匹配.然而这些方法均依赖于精确提取的二值化图像, 在环境光照度变化, 长时间观测时小鼠对背景环境不断产生影响的条件下, 二值化预处理过程中均不可避免地引入噪声干扰, 进而降低了上述体态识别算法的准确率[15], 限制了这些方法的应用环境.
事实上, 实验小鼠运动参数提取的主要困难源于其非刚性变形的特征.在鼠笼等被观测视野较小的条件下, 这一特点尤为显著.如果再考虑视频信号的干扰, 则采用肢体的不变特征用于连续帧间的匹配存在较大的困难[16].为解决这一困境, 本文引入模板匹配和粒子滤波的基本识别框架, 提出实验小鼠的几何体部件模型及视频观测与模型匹配方法, 借助粒子滤波器完成噪声干扰条件下的小鼠运动轨迹参数及体态参数的同时提取.
1. 目标跟踪方法框架
实验小鼠运动参数的提取问题本质上可以归结为一个连续视频中的目标跟踪问题.而目标跟踪在安全监控、气象分析、医疗诊断、军事制导等领域已经具备大量研究成果[17], 其核心问题是在后续图像中定位和搜索待跟踪目标的先验特征.目标跟踪的关键因素包括对象的表示方法、搜索算法、模型更新等, 现已形成了主动轮廓跟踪、基于特征的跟踪、基于模型的跟踪以及基于判别式分类的跟踪等几类方法[18].本文所采用的模板匹配和粒子滤波方法是视频跟踪的一种重要方法, 其基本原理如图 1所示.
图 1 用于小鼠目标跟踪的模板匹配及粒子滤波方法Fig. 1 Model matching and particle filtering method applied to the object tracking of a laboratory mouse图 1中的观测结果是待处理图像的二值化结果, 视野内标记为黑色的像素是小鼠目标及干扰像素. 图 1中的粒子集合是小鼠诸多可能状态所对应的理想观测结果.将实际观测结果与之进行对比匹配, 并从中挑选出匹配度最佳的运动状态作为观测结果.经过粒子滤波重采样后, 生成了一批新的状态粒子以备下一帧图像的匹配处理.
模板匹配方法避免了直接从精确二值化图像中拟合小鼠特征参数的困境, 其优势是降低了图像预处理的难度, 并允许一定程度的图像噪声干扰存在, 因而对系统的前期参数调整和标定要求较低.同时, 模板匹配方法不对二值化实施过程做出过多假设, 无论是阈值划分、背景减除或是多种策略的综合方法, 只要能够尽可能凸显小鼠目标和活动背景的差异, 即可在统一框架下实现运动参数的提取, 因而扩展了方法对光照、背景等条件变化的适应性.
基于模板匹配和粒子滤波的关键是建立小鼠外形的观测模型以及模拟小鼠活动的运动模型, 以前者作为粒子滤波器的观测模型, 以后者作为状态转移模型, 从而建立起如图 1所示的完整目标跟踪算法.
2. 实验小鼠的物理模型
小鼠的物理模型与观测视角相关.一般有顶置摄像机视角[7-11]、侧置摄像机视角[3, 14]、斜置摄像机视角[5]及各种混合摄像机视角等.其中顶置摄像机能够观测小鼠的运动轨迹和移动速度, 也是观测体态的较优视角, 因此物理模型主要针对顶置摄像机观测情形.
2.1 实验小鼠模型
小鼠的模型如图 2所示:将小鼠分解为头部和臀部, 从而得到两个椭圆形几何体构成的部件模型.两个椭圆形的几何尺寸固定不变, 但两者的相对位置发生变化. 图 2中, $Oxy$ 是模型坐标系, 固连在臀部椭圆形的中心.头部椭圆的中心点表示为 $c_h$ , 头部相对臀部的运动轨迹表示为 $c_h$ 沿 $C$ 点为圆心的一段圆弧, $C$ 点仅位于 $Oxy$ 的 $y$ 轴上, 坐标为 $( {0,y_c } )$ , 定义体态的曲率为
$\begin{eqnarray} \rho = \frac{1}{y_{c}},\quad {y_{c}} \ne 0 \end{eqnarray}$
(1) 与一般圆弧曲率定义不同, $\rho$ 带有符号, 当 $\rho<0$ , 圆心 $C$ 位于 $x$ 轴下方, 头部以顺时针方向偏转, 当 $\rho>0$ , 圆心 $C$ 位于 $x$ 轴上方, 头部向逆时针方向偏转, 当 $\rho=0$ 时, 圆弧退化为直线, 此时头部相对臀部仅有伸缩运动, 若设圆弧长 $e>0$ , 则圆弧对应的圆心角可表示为 $\phi =e\rho$ .
头部椭圆的半长轴、半短轴分别是 $(a_{h},b_{h})$ , 臀部椭圆的半长轴、半短轴分别是 $(a_{t},b_{t})$ . $\rho =0$ , $e=0$ 时, 模型处于零点, 此时 $c_{h}$ 的坐标为 $\left({f,0} \right)$ . $a_{h},b_{h} ,a_{t},b_{t},f$ 与被跟踪小鼠的体型有关, 在视频跟踪前通过手动设定, 并在跟踪过程中保持不变.因此, 小鼠的体态可由体态曲率 $\rho$ 和体态伸长 $e$ 两个参数描述.
上述模型本质上是一种变形体的近似部件模型[18-19].依照这一思路建立的其他模型包括椭圆包围盒模型[13]及多几何元素模型[12]等.相比于这些已有模型, 本文提出新模型主要考虑了几个因素的综合与平衡, 具体包括外形相似性、参数可视化和计算高效性等三个方面.
1) 新模型将小鼠简化为头部、臀部两个部分, 能够更为细致地描绘蜷缩和弯曲等体态.而单个椭圆包围盒模型则以其长短轴作为描述参数, 仅能有效描述体长参数, 而无法描述体态曲率.另外, 从生物学角度分析实验小鼠的骨骼特征, 其体态的变化对应了脊柱的状态, 而模型中圆弧实现了脊柱状态在顶置摄像机视角方向的近似变化.尽管文献[13]中的模型从机械连接及物理学规律应用的角度较为容易, 但缺乏这一相似性.
2) 新模型可将小鼠体态的关键参数归纳为体态伸长 $e$ 和体态曲率 $\rho$ .同时将臀部(或模型固连坐标系 $Oxy$ )的位姿及移动速率描绘为小鼠整体的运动状态, 这些参数完整展现了先前研究过程中小鼠行为学分析的几个关键参数[11, 13], 从而在参数可视化方面具有优势.
3) 在给定模型状态的条件下, 只需要针对二值化图像中的前景像素点执行两个椭圆形的包含性测试即可完成匹配检验和评估, 从而具有更低的计算量.同时, 两个椭圆间仅存在圆弧的相对运动, 因而两者相对位置计算较为容易.这些特征对后续粒子滤波跟踪过程中大量的匹配计算性能尤为重要.
当然, 模型也存在几个尚未考虑问题:其一臀部椭圆模型尺寸固定无法反应蜷缩时臀部形状的真实变化, 其二是两个椭圆相交处并非光滑过渡, 与实际小鼠的外形略有差异.但是这些参数对小鼠行为学研究的意义较小, 如果引入更多的状态参数或局部部件将大幅度增加滤波提取的难度.因而这些因素就作为算法的系统性偏差或者未来用于在线模型修正的因素而暂时不再考虑.
2.2 小鼠运动状态方程
通过建立实验小鼠的运动状态方程以实现其运动状态的描述, 两个关键问题是状态向量中各个分量的选择以及状态方程的建立.
针对第2.1节中的小鼠部件模型，其完整运动状态应包含两部分:体态变化及整体位姿变化.前者由体态伸长 $e$ 和体态曲率 $\rho$ 两个变量所描述, 后者则包括平面位姿变量 $\left({x,y,\theta}\right)$ 表示其在全局坐标系下的坐标及姿态角, 从而构成基本状态向量.基于小鼠运动不受外界主观控制的特点, 一种运动建模方法是假定各分量服从独立随机变化规律的随机运动模型, 即在前一个状态基础上增加一个随机偏移从而推算出新状态, 结果如式(2)所示.
$\begin{eqnarray} {\pmb X}_{k+1} ={\pmb X}_k +{\boldsymbol \delta } \end{eqnarray}$
(2) 然而这种方法没有充分利用实验小鼠的运动习性特征, 小鼠的运动状态变为纯粹猜测, 最终影响了精确估计条件下的滤波算法的粒子规模.同时小鼠的移动速度仍然需要对小鼠位置变量 $\left({x,y}\right)$ 做差分运算, 丧失了粒子滤波充分利用先验观测实现精确估计的优势.
为解决上述问题, 在状态向量中引入一个移动速率变量 $v$ , 并认为 $v$ 的方向与头部椭圆模型的长轴方向一致.如图 3所示, 在移动速率 $v$ 及体态曲率 $\rho$ 恒定不变的条件下, 小鼠将以匀速圆周规律运动.其圆心与图 2中的 $C$ 点重合, 运动轨迹的曲率则与式(1)中的 $\rho$ 相同.
图 3 小鼠以恒定的体态曲率 $\rho$ 及移动速率 $v$ 的匀速圆周运动Fig. 3 Uniform circular motion of a mouse with constant curvature $\rho $ and constant moving velocity $v$由于小鼠体态的变化难以预测, 因此体态伸长 $e$ 和体态曲率 $\rho$ 仍服从随机运动规律.综上所述, 新的状态向量可表示为
$\begin{eqnarray*} { {\pmb X}}=\left( {x,y,\theta ,v,e,\rho } \right)^{\rm T} \end{eqnarray*}$
而状态方程则如式(3)所示.
$\begin{eqnarray} \left( {\begin{array}{*{20}{c}} {{x_{k + 1}}}\\ {{y_{k + 1}}}\\ {{\theta _{k + 1}}}\\ {{v_{k + 1}}}\\ {{e_{k + 1}}}\\ {{\rho _{k + 1}}}\\ \end{array} } \right) = \left( {\begin{array}{*{20}{c}} {{f_x}({\theta _k},{v_k},{\rho _k},{x_k},{y_k})}\\ {{f_y}({\theta _k},{v_k},{\rho _k},{x_k},{y_k})}\\ {{\theta _k} + {v_k}{\rho _k}\Delta t}\\ {{v_k}}\\ {{e_k}}\\ {{\rho _k}}\\ \end{array} } \right){\text{ + }}{\boldsymbol \delta} \end{eqnarray}$
(3) 式中, $\Delta t$ 是视频采集周期, 函数 $f_x$ 与 $f_y$ 分别是小鼠位置 $x$ 和 $y$ 的变化规律, 基于圆周运动假设有:
$\begin{align} &\left( {\begin{array}{*{20}{c}} {{f_x}({\theta _k},{v_k},{\rho _k},{x_k},{y_k})} \\ {{f_y}({\theta _k},{v_k},{\rho _k},{x_k},{y_k})} \\ \end{array} } \right) = \nonumber\\ &\qquad \left( {\begin{array}{*{20}{c}} {\cos {\theta _k}} & { - \sin {\theta _k}} \\ {\sin {\theta _k}} & {\cos {\theta _k}} \\ \end{array} } \right)\left( {\begin{array}{*{20}{c}} {\sin {\omega _k}\Delta\frac{ t}{\rho _k}} \\ \frac{(1 - \cos {\omega _k}\Delta t)}{{\rho _k}} \\ \end{array} } \right) + \nonumber\\ &\qquad\left( {\begin{array}{*{20}{c}} {{x_k}} \\ {{y_k}} \\ \end{array} } \right) \end{align}$
(4) 其中, $\omega _k $ 是圆周运动的角速度, 且有.式(3)中， ${\rm {\boldsymbol \delta }}$ 是六维随机向量, 描述了小鼠运动的随机特性.
此外, 实际应用状态方程(3)时, 还需考虑小鼠体态的物理限制、移动速度的限制以及视频观测的视野范围等, 针对式(3)所得的状态向量 ${ {\pmb X}}$ 的每一个分量 $q$ 限幅:
$\begin{eqnarray} \begin{array}{*{20}{lll}} {{q_i} = \left\{ {\begin{array}{*{20}{lll}} {q_i^{\max }},& {q_i^{\max } < {q_i}}\\ {{q_i}},& {q_i^{\min } \leq {q_i} \leq q_i^{\max }},\\ {q_i^{\min }},& {{q_i} < q_i^{\min }} \end{array} } \right.} ~&~ {i = 1,2,\cdots,7} \\ \end{array} \end{eqnarray}$
(5) 式中, $q_1,q_2 ,\cdots,q_7 $ 分别对应了的各个分量 $x,y,\theta,v,e,\rho $ , 各变量物理限制依据经验预先设定.
3. 实验小鼠的粒子滤波跟踪
粒子滤波是Bayes滤波器的一种形式, Bayes滤波器的目标是给出状态变量在一次新观测条件下的后验概率分布.粒子滤波器通过蒙特卡洛方法, 以一定数量的采样粒子来近似表示这一概率分布, 其特点是能够适用于非线性系统, 并易于理解和实现.粒子滤波算法已具有坚实的理论基础和成熟的技术实现, 细节可参照文献[20].本文采用基本粒子滤波算法, 着重讨论实验小鼠跟踪这一特定问题下应考虑的因素, 至于针对粒子匮乏、粒子抢夺等缺陷的各种算法改进和提升可参考文献[21].
3.1 基本粒子滤波算法
基本粒子滤波算法列表如表 1所示.
1: Function Particle_filter $\chi _t $ , $z_{t+1}$ 2: 3: For $m = 1:M$ 4: 以采样获得 5: 计算权重 6: 7: End for 8: 依据 $\bar {\chi }_{t+1} $ 重采样生成 $\chi _{t+1} $ 9: Return $\chi _{t+1} $ 10: End function 算法中 $M$ 是预设的粒子总数, $\chi _t $ 是上一个周期的 $M$ 个粒子集合, $z_{t+1} $ 是本周期的观测结果, 对于小鼠跟踪这一特定问题, 是图像的二值化图像结果.第4行的采样过程即依据小鼠运动模型式(3)从前一个状态产生新的状态 ${\rm {\pmb X}}_{t+1}^{[m]} $ .第5行计算对应的理论观测结果和实际观测结果 $z_{t+1} $ 两者的匹配度, 从而求得状态粒子的重要性权重 $w_{t+1}^{[m]} $ , 是所有粒子状态及其权重二元组的集合.第8行则是粒子滤波器的重采样过程, 由于重采样算法较为成熟, 且与小鼠跟踪这一特定问题较为独立, 这里不再详细列出, 具体可参照文献[20-21].
由此可见, 为实现上述粒子滤波算法, 需要为式(3)中的随机向量定义合理的概率分布规律.同时需要依据实际观测 $z_{t+1} $ 求出任意粒子状态的重要性权重 $w_{t+1}^{[m]} $ .
3.2 模板匹配及重要性权重计算方法
设摄像机的采集分辨率为 $W\times H$ 像素, 则算法列表 1中第5行中的 $z_{t+1} $ 由个布尔函数 $q(u,v)$ 组成, 其中, $(u,v)$ 是像素坐标, $q(u,v)=\text{true}$ 表明该像素隶属于小鼠目标; 反之, 则不属于小鼠目标, 记所有属于小鼠的像素坐标集合为 $\Sigma _O $ , 即:
$\begin{eqnarray} \Sigma _O =\{(u,v)\vert q(u,v)=\text{true}\} \end{eqnarray}$
(6) 另一方面, 若将所对应的小鼠几何体模型绘制在 $W\times H$ 的空白二值图像中, 并计 $q(u,v)=\text{true}$ 为几何模型的内部像素, $q(u,v)=\text{false}$ 为几何模型的外部像素.则可以获得与式(6)类似的小鼠部件模型的像素坐标集合, 若记 $\Sigma _{O} $ 与中完全匹配的像素点集合为 $\Sigma _{M}^{[m]} $ , 则可表示为
$\begin{eqnarray} \Sigma _{M}^{[m]} =\{p\vert p\in \Sigma _{O} ,p\in \Sigma _{P}^{[m]} \} \end{eqnarray}$
(7) 为描述实际观测与理论模型之间的差异, 定义一种差异度函数为
$\begin{eqnarray} \begin{split} (r^{[m]})^2 =~& \frac{1}{2}\left[{\frac{n(\Sigma _{P}^{[m]} )-n(\Sigma _{M}^{[m]} )}{n(\Sigma _{M}^{[m]} )}} \right]^2+ \\ & \frac{1}{2}\left[{\frac{n(\Sigma _{O} )-n(\Sigma _{M}^{[m]} )}{n(\Sigma _{M}^{[m]} )}} \right]^2 \\ \end{split} \end{eqnarray}$
(8) 其中, $n(\cdot )$ 为集合中元素的个数.当时, 由式(8)计算出 $r^{[m]}=0$ , 表示实际观测与该状态粒子的差异度为0, 两者完全匹配; 当时, 由式(8)可知, 表示实际观测与部件模型不具备任何匹配像素, 两者完全不同.而当 $\Sigma _{M}^{[m]} $ 的元素数目从0逐步增加至时, 逐步减小至0, 从而表明两者差异度逐步降低直至完全匹配.式(8)连续可微且对匹配像素数 $n(\Sigma_{M}^{[m]})$ 单调递减, 因而可用于生成粒子的重要性权重系数.
将 $r^{[m]}$ , $1\le m\le M$ 视作一个随机变量 $r$ 的一系列采样, 并假设 $r$ 服从均值为0方差为 $\sigma_0^2 $ 的高斯分布, 则 $r$ 在 $r^{[m]}$ 处对应的概率密度为
$\begin{eqnarray} w_{t+1}^{[m]} =N\left( {r^{[m]},0,\sigma _0^2 } \right) \end{eqnarray}$
(9) 其中, $N\left( {r,0,\sigma _0^2 } \right)$ 为高斯分布的概率密度函数, 由于该概率密度函数相对于 $r=0$ 轴对称, 因此有:
$\begin{eqnarray} w_{t+1}^{[m]} =N\left( {\left| {r^{[m]}} \right|,0,\sigma _0^2 } \right) \end{eqnarray}$
(10) 结合式(8)和式(10)得到第 $m$ 个粒子的差异度 $r^{[m]}$ 的概率密度, 作为第 $m$ 个样值粒子所对应的相对权重, 并用于后续重采样过程, 从而得到算法列表 1第5行的粒子权重计算方法.
3.3 运动模型中的随机参数
表 1中第4行的计算依据是式(3)所示的实验小鼠运动状态方程.粒子状态的推移与式中六维随机向量的分布规律有关, 本文假定其服从均值为 ${\rm {\bf 0}}$ , 协方差矩阵为 ${\boldsymbol \Sigma} $ 的六维联合Guass分布.然而相应的随机实验存在较大困难, 无法确定的具体数值.为简化计算难度, 进一步忽略了随机向量的各个分量的相关性, 从而将六维随机向量简化为服从相互独立的Guass分布, 其协方差矩阵 ${\boldsymbol \Sigma} $ 退化为一个对角阵, 对角线元素分别是每一个随机变量的方差, 即:
$\begin{eqnarray} \boldsymbol \Sigma = \text{diag}\{\sigma _x^2,\sigma _y^2,\sigma _\theta ^2,\sigma _v^2,\sigma _e^2,\sigma _\rho ^2\} \end{eqnarray}$
(11) 3.4 粒子的密度估计
粒子滤波器以一系列粒子状态 $\chi _{t+1} $ 模拟小鼠的状态的概率分布, 然而最终需要由 $\chi _{t+1} $ 给出状态的最佳估计作为最终求解结果.该过程是一个密度估计过程, 具体方法包括均值估计、 $K$ -means聚类、直方图法以及核密度估计法等[20].
诸多方法中最易实施、计算代价最低的是均值估计, 即假设粒子分布满足高斯分布, 计算各粒子对应状态的期望和方差, 并将期望作为估计结果.然而该方法无法处理小鼠在透明鼠笼上的镜像所导致的多峰分布状况, 均值运算将混合小鼠目标与其镜像状态, 从而可能产生较大的跟踪误差.核密度估计法能够拟合出一批采样粒子所对应的精确分布规律, 但是多元核密度估计算法较为复杂, 且核密度估计后还需求解非线性函数的局部极值才能获得最终结果, 计算效率较低.直方图法则难以应对六维状态空间的有效划分问题, 难以在有限的存储空间下获得优良的估计结果.
相比之下, $K$ -means聚类法[22]依据状态向量间的欧氏距离对采样粒子进行分类, 可以满足多峰值条件下的提取要求.同时, 状态向量的欧氏距离计算能够满足计算效率要求, 因而在多个多方面具有优势. $K$ -means算法在实际应用过程中需要预先指定峰值的最多个数, 其来源于对小鼠状态分布的先验知识.例如, 粒子状态的峰值主要出现于小鼠首尾对称的两种状况下, 则可预先指定聚类类别 $k=2$ , 若再考虑实验鼠笼产生小鼠镜像的状况, 则可设定 $k=4$ . $K$ -means聚类后, 从最终峰值中, 选择权重之和最大的一个作为最终的识别结果.
4. 实验结果及讨论
4.1 实验装置及方法
为验证小鼠运动参数提取算法的有效性, 设计并实施了一种小动物行为学实验装置, 系统原理如图 4所示.主体结构是一个铁制白色不透光箱体, 内部可放置若干个上端开口的透明动物笼.箱体四壁分布了多个LED照明灯管.照明自动控制系统实现了三种照明模式:自然强度的白光、强照度蓝光及黑暗条件下的红外光.三种模式可以以任意次序, 任意时长组合并循环切换.箱体顶部配置一台CCD彩色摄像机, 输出的信号经视频采集系统转换为数字图像并传送到图像处理工作站, 经工作站在线处理并保存为AVI格式视频文件, 视频采集帧速率为24帧/s.图形处理工作站采用Intel Core i7-4700MQ处理器及8 GB内存.
图 4 带有视频跟踪的小鼠行为学实验装置Fig. 4 Experimental device for behavior analysis of mice with video tracking行为学实验要求在多种光照模式下观测并提取小鼠的运动轨迹、运动速度和体态等参数, 实验时间最长可达48 h.采用白色小鼠在三种光照条件下的样本图片如图 5所示.
4.2 小鼠跟踪实验及结果
为了适应不同的光照条件, 避免频繁的颜色域值标定, 采用背景减除法提取小鼠的二值化图像.分别在不同的光照条件下采用混合Guass模型对一段初始图像序列执行背景学习后, 启动背景减除过程, 同时在背景减除处理期间, 每间隔若干帧执行一次背景学习, 从而不断适应环境所发生的微小变化.对背景减除所生成的二值化图像进行了形态学运算和轮廓求解, 从而滤除了明显的识别噪声.上述过程采用OpenCV图像处理库实现.蓝色强光照射条件下的一帧处理实例如图 6所示.由图 6可见, 最终处理结果除了目标小鼠的像素外, 还遗留了透明鼠笼所产生的镜像等干扰尚未滤除.
图 6 背景减除及后续的图像预处理过程示例Fig. 6 An example of the background subtraction and the followed image preprocessing粒子滤波算法基于Bayes滤波库OROCOS [23]实现, 所形成的算法流程图如图 7所示.
图 7 图像预处理与粒子滤波的小鼠运动参数提取流程图Fig. 7 Flow chart of the parameters for behavior analysis based on image preprocessing and particle filter图 7所示的过程中, 粒子滤波器的粒子总数设定为500, 当无效粒子权重达到粒子总数的1/4启动重采样.实验过程中, 式(3)所对应的物理量如表 2所示, 此外观测模型中的模型差异度 $r$ 及其概率分布参数也列在表 2中.
表 2 实验小鼠状态变量及模型差异度 $r$ 的正态分布参数Table 2 Parameters of the normalized distribution for the state variables and model-observation difference $r$ of a laboratory mouse变量 单位 均值 $\mu $ 标准差 $\sigma $ 位置 $x$ 像素 $0.00$ $3.00 $ 位置 $y$ 像素 $0.00$ $3.00 $ 姿态角 $\theta $ rad $0.00$ $0.09 $ 速度 $v$ 像素 $\cdot \text{s}^{-1}$ $0.00$ $25.00 $ 体态伸长 $e$ 像素 $0.00$ $3.00 $ 体态曲率 $\rho $ $\text{rad}\cdot$ 像素 $^{-1}$ $0.00$ $0.01 $ 模型差异度 $r$ 1 $0.00$ $1.00 $ 如图 8所示是在蓝色强光照射条件下截取的一段小鼠视频跟踪结果.在透明鼠笼镜面反射干扰的条件下, 能够有效地跟随小鼠的运动轨迹, 除较为精确的位置 $x,y$ 外, 还包括其姿态角 $\theta $ , 且能够有效地区分小鼠的头部和尾部.
图 8 蓝光照射条件下的一段小鼠目标跟踪结果Fig. 8 One piece of the object tracking result of a mouse in the blue illuminating condition以蓝光-白光-黑暗三种照射条件各持续60 s为标准实验, 每一种光照条件下的前20 s用于背景学习, 后40 s启动背景减除及粒子滤波算法.将三段的提取结果进行连接, 小鼠的运动轨迹和姿态角如图 9所示, 可见其轨迹连续、姿态角和轨迹具有良好的一致性, 未出现因为识别错误而导致的小鼠位置跳跃及轨迹不连续现象.因而证实本文的方法能够达到优良整体轨迹跟踪性能.
图 9 模板匹配及粒子滤波方法在三种不同光照条件下所提取的实验小鼠整体运动轨迹提取结果Fig. 9 The extracting result of the moving trajectory of the laboratory mouse in three different illuminating conditions using model matching and particle filtering如图 10所示是算法所提取的运动参数, 包括小鼠移动速度图 $v$ 、体态伸长 $e$ 和体态曲率 $\rho $ 变化图.其能够较为精确地显示出实验小鼠的运动速度和体态的变化. 图 10中, 在实验小鼠具有较高移动速度时, 均对应了其体态伸长幅度的增加, 因而符合小鼠的运动习性.其中65 s $\sim $ 80 s时间段内, 小鼠的移动速度接近于零, 同时体态伸长也缩减到最小值表明小鼠处于静止蹲窝的状态, 与其他时间段内小鼠的快速移动和探索行为对比明显, 因而所提取的小鼠体态及其他运动参数可用于小鼠的后续行为学分析.
图 10 模板匹配及粒子滤波算法所提取的移动速度、体长变化及体态曲率结果Fig. 10 The extracting result of the motion speed, body length variation and body curvature using model matching and particle filtering上述的三种光照切换实验中, 排除前期背景学习实验段的每一独立帧的处理时间统计直方图如图 11所示, 可见该时间主要分布在30 ms以下, 因而能够达到实验装置24帧/s的采集速率下的实时在线提取需求.
图 11 单周期运动参数提取算法计算耗时直方图Fig. 11 Histogram of the computing time consumption of the motion parameter extraction algorithm for a single cycle4.3 与单帧图像差分跟踪方法的对比实验
在本文所述的模板匹配和粒子滤波算法外, 还采用面向单帧数字图像处理的方法求解二值化图像中前景像素的矩形包围盒, 并通过帧间差分的方法进行了实验小鼠轨迹跟踪以及移动速度提取的实验.如图 12所示是基于相同的背景减除方法及形态学滤波及轮廓计算, 针对如图 9所示相同的实验视频的轨迹提取结果, 图 13是利用差分法计算所得的移动速度和体长变化参数.
图 12 采用帧间差分法在三种光照条件下所提取的目标小鼠整体运动轨迹提取结果Fig. 12 The extracting result of the moving trajectory of the object mouse in three different illuminating conditions using frame differencing method
图 13 帧间差分法所提取的移动速度、体长变化及体态曲率结果Fig. 13 The extracting result of the motion speed, body length variation, and body curvature using frame differencing method图 12中的轨迹存在较大的跳跃, 并且存在一系列误识别结果, 这些错误的识别结果来源于将实验小鼠的镜像误当做小鼠本身.同时图 13对小鼠移动速度和体长参数的估计被大量噪声所干扰, 难以统计出活动频率的特征, 进而很难应用于后续行为参数的分析.
4.4 针对移动速度 $\pmb v$ 的讨论
如第2.2节中的状态向量中引入小鼠移动速率 $v$ 的策略除了能够实现该变量的直接估计外, 还能够在目标跟踪时尽快纠正首尾倒置的错误.
实验表明, 当图 2所示的模型中臀部椭圆参数与头部椭圆参数相近或者图像二值化过程中存在较大识别误差的条件下, 易导致小鼠首尾倒置的识别结果.此时若按照式(2)的随机运动规律, 则可能在较长时间内无法纠正这一错误.
考虑到小鼠更习惯于前进运动的习性, 在引入移动速率 $v$ 且在式(5)的限幅函数中设定小鼠的最高前进移动速率远大于最高后退移动速率的条件下, 只有那些首尾状态与移动速率方向一致的粒子才能够达到更好的匹配度, 而首尾状态与移动速度方向相反的粒子则迅速遭到淘汰, 从而能够在后续几帧内快速消除首尾的识别错误, 提高了小鼠运动的跟踪准确度.
4.5 粒子滤波器的必要性讨论
本文所提出的方法能够有效滤除透明鼠笼产生的虚影, 并得出较为精确的参数提取结果得益于粒子滤波器的几项特征.粒子滤波基于概率随机思想, 求解结果是小鼠状态的概率分布而非单一取值, 当概率分布存在多个峰值时, 可以通过小鼠的活动范围对结果进行进一步过滤, 从而排除了错误的识别结果.同时, 基于Markov过程的方法能够有效利用多次观测结果, 逐步改正识别错误和观测误差, 从而避免错误估计的传播.粒子滤波器本质上建立了视频图片的前后联系, 协助滤除了图像二值化的噪声, 从而在较大干扰的条件下, 依然能够实现小鼠宏观运动和微观体态信息的同时精确提取.
5. 结论
本文将实验小鼠的运动参数提取归结为一个视频目标跟踪问题, 并采用模板匹配和粒子滤波思想对该问题进行求解.提出一种实验小鼠的简化几何体模型, 在完整描绘其运动特性的前提下, 充分考虑了模型与顶置视角实验小鼠的形体相似性, 并实现了高效的匹配运算性能.所提出的运动模型中引入了小鼠移动速率参数, 在实现了这一参数的精确估计之外, 有效解决了首尾识别混淆的错误, 提升了参数提取的精确度.粒子滤波方法有效建立起了视频图像的前后联系, 解决了虚影等一系列干扰问题.与逐帧识别差分法的对比实验证实了模板匹配与粒子滤波器能够实现强噪声干扰条件下小鼠宏观运动和微观体态的同时提取, 从而能够为后续行为学分析提供可靠依据.
未来研究方向包括:依据识别结果更新小鼠模型长短轴参数的在线自适应调整; 基于粒子滤波器具有融合多个视频采集信息的能力, 实现更为精确的运动参数乃至三维运动参数的提取, 可能的解决方案如建立一种合理的小鼠三维观测模型; 此外如果考虑多只小鼠的交互行为实验, 则需要对粒子滤波器做进一步的改进, 以解决高维度下的运算负荷问题和粒子抢夺等问题.
-
图 1 用于小鼠目标跟踪的模板匹配及粒子滤波方法
Fig. 1 Model matching and particle filtering method applied to the object tracking of a laboratory mouse
图 3 小鼠以恒定的体态曲率 $\rho$ 及移动速率 $v$ 的匀速圆周运动
Fig. 3 Uniform circular motion of a mouse with constant curvature $\rho $ and constant moving velocity $v$
图 4 带有视频跟踪的小鼠行为学实验装置
Fig. 4 Experimental device for behavior analysis of mice with video tracking
图 6 背景减除及后续的图像预处理过程示例
Fig. 6 An example of the background subtraction and the followed image preprocessing
图 7 图像预处理与粒子滤波的小鼠运动参数提取流程图
Fig. 7 Flow chart of the parameters for behavior analysis based on image preprocessing and particle filter
图 8 蓝光照射条件下的一段小鼠目标跟踪结果
Fig. 8 One piece of the object tracking result of a mouse in the blue illuminating condition
图 9 模板匹配及粒子滤波方法在三种不同光照条件下所提取的实验小鼠整体运动轨迹提取结果
Fig. 9 The extracting result of the moving trajectory of the laboratory mouse in three different illuminating conditions using model matching and particle filtering
图 10 模板匹配及粒子滤波算法所提取的移动速度、体长变化及体态曲率结果
Fig. 10 The extracting result of the motion speed, body length variation and body curvature using model matching and particle filtering
图 11 单周期运动参数提取算法计算耗时直方图
Fig. 11 Histogram of the computing time consumption of the motion parameter extraction algorithm for a single cycle
图 12 采用帧间差分法在三种光照条件下所提取的目标小鼠整体运动轨迹提取结果
Fig. 12 The extracting result of the moving trajectory of the object mouse in three different illuminating conditions using frame differencing method
图 13 帧间差分法所提取的移动速度、体长变化及体态曲率结果
Fig. 13 The extracting result of the motion speed, body length variation, and body curvature using frame differencing method
1: Function Particle_filter $\chi _t $ , $z_{t+1}$ 2: $\bar {\chi }_{t+1} \leftarrow \emptyset,\chi _{t+1} \leftarrow \emptyset $ 3: For $m = 1:M$ 4: 以 ${\rm {\pmb X}}_{t+1}^{[m]} \sim p({\rm {\pmb X}}_{t+1} \vert {\rm {\pmb X}}_t )$ 采样获得 ${\rm {\pmb X}}_{t+1}^{[m]} $ 5: 计算权重 $w_{t+1}^{[m]} \mbox{=}p(z_{t+1} \vert {\rm {\pmb X}}_{t+1}^{[m]} )$ 6: $\bar {\chi }_{t+1} =\bar {\chi }_{t+1} +\left\langle {{\rm {\pmb X}}_{t+1}^{[m]} ,w_{t+1}^{[m]} } \right\rangle $ 7: End for 8: 依据 $\bar {\chi }_{t+1} $ 重采样生成 $\chi _{t+1} $ 9: Return $\chi _{t+1} $ 10: End function 
下载: 导出CSV
表 2 实验小鼠状态变量及模型差异度 $r$ 的正态分布参数
Table 2 Parameters of the normalized distribution for the state variables and model-observation difference $r$ of a laboratory mouse
变量 单位 均值 $\mu $ 标准差 $\sigma $ 位置 $x$ 像素 $0.00$ $3.00 $ 位置 $y$ 像素 $0.00$ $3.00 $ 姿态角 $\theta $ rad $0.00$ $0.09 $ 速度 $v$ 像素 $\cdot \text{s}^{-1}$ $0.00$ $25.00 $ 体态伸长 $e$ 像素 $0.00$ $3.00 $ 体态曲率 $\rho $ $\text{rad}\cdot$ 像素 $^{-1}$ $0.00$ $0.01 $ 模型差异度 $r$ 1 $0.00$ $1.00 $
下载: 导出CSV
-
[1] 尚玉昌.动物行为学.北京:北京大学出版社, 2005.Shang Yu-Chang. Animal Ethology. Beijing:Peking University Press, 2005. [2] Dell A I, Bender J A, Branson K, Couzin I D, de Polavieja G G, Noldus L P J J, Pérez-Escudero A, Perona P, Straw A D, Wikelski M, Brose U. Automated image-based tracking and its application in ecology. Trends in Ecology & Evolution, 2014, 29(7):417-428 https://www.researchgate.net/publication/263087352_Automated_image-based_tracking_and_its_application_in_ecology [3] Jhuang H, Garrote E, Yu X L, Khilnani V, Poggio T, Steele A D, Serre T. Automated home-cage behavioural phenotyping of mice. Nature Communications, 2010, 1(6):Article No.68 https://dspace.mit.edu/openaccess-disseminate/1721.1/64492 [4] Ishii H, Ogura M, Kurisu S, Komura A, Takanishi A, Iida N, Kimura H. Development of autonomous experimental setup for behavior analysis of rats. In:Proceedings of the 2007 IEEE/RSJ International Conference on Intelligent Robots and Systems. San Diego, USA:IEEE, 2007. 4152-4157 [5] Salem G, Dennis J, Krynitsky J, Garmendia-Cedillos M, Swaroop K, Malley J D, Pajevic S, Abuhatzira L, Bustin M, Gillet J P, Gottesman M M, Mitchell J B, Pohida T J. SCORHE:a novel and practical approach to video monitoring of laboratory mice housed in vivarium cage racks. Behavior Research Methods, 2014, 47(1):235-250 doi: 10.3758/s13428-014-0451-5.pdf [6] Noldus L P J J, Spink A J, Tegelenbosch R A J. EthoVision:a versatile video tracking system for automation of behavioral experiments. Behavior Research Methods, Instruments, & Computers, 2001, 33(3):398-414 doi: 10.3758/BF03195394 [7] Rantalainen T, Silvennoinen M, Kainulainen H, Sievänen H. Vertical ground reaction force measurements and video measurements provide comparable estimates of distance moved by mice during artificial light and dark periods. Journal of Neuroscience Methods, 2011, 197(1):104-108 doi: 10.1016/j.jneumeth.2011.02.010 [8] Dielenberg R A, Halasz P, Day T A. A method for tracking rats in a complex and completely dark environment using computerized video analysis. Journal of Neuroscience Methods, 2006, 158(2):279-286 doi: 10.1016/j.jneumeth.2006.05.024 [9] Zurn J B, Jiang X H, Motai Y. Video-based tracking and incremental learning applied to rodent behavioral activity under near-infrared illumination. IEEE Transactions on Instrumentation and Measurement, 2008, 56(6):2804-2813 https://www.researchgate.net/publication/3094058_Video-Based_Tracking_and_Incremental_Learning_Applied_to_Rodent_Behavioral_Activity_Under_Near-Infrared_Illumination [10] Pistori H, Odakura V V V A, Monteiro J B O, Gonçalves W N, Roel A R, de Andrade Silva J. Mice and larvae tracking using a particle filter with an auto-adjustable observation model. Pattern Recognition Letters, 2010, 31(4):337-346 doi: 10.1016/j.patrec.2009.05.015 [11] 张敏, 张恒义, 郑筱祥.基于轮廓曲率和谱系聚类的大鼠体态自动识别.浙江大学学报(工学版), 2006, 40(3):524-527, 532 https://www.wenkuxiazai.com/doc/86e34c07227916888486d75c.htmlZhang Min, Zhang Heng-Yi, Zheng Xiao-Xiang. Automatic recognition of rat's postures based on contour curvature and hierarchical clustering analysis. Journal of Zhejiang University (Engineering Science), 2006, 40(3):524-527, 532 https://www.wenkuxiazai.com/doc/86e34c07227916888486d75c.html [12] de Chaumont F, Coura R D S, Serreau P, Cressant A, Chabout J, Granon S, Olivo-Marin J C. Computerized video analysis of social interactions in mice. Nature Methods, 2012, 9(4):410-417 doi: 10.1038/nmeth.1924 [13] Shi Q, Miyagishima S, Fumino S, Konno S, Ishii H, Takanishi A. Development of a cognition system for analyzing rat's behaviors. In:Proceedings of the 2010 IEEE International Conference on Robotics and Biomimetics (ROBIO). Tianjin, China:IEEE, 2010. 1399-1404 [14] Branson K, Belongie S. Tracking multiple mouse contours (without too many samples). In:Proceedings of the 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. San Diego, CA, USA:IEEE, 2005:1039-1046 [15] Farah R, Langlois J M P, Bilodeau G. Catching a rat by its edglets. IEEE Transactions on Image Processing, 2013, 22(2):668-678 doi: 10.1109/TIP.2012.2221726 [16] de Andrade Silva J, Gonçalves W, Machado B, Pistori H, de Souza A S, de Souza K P. Comparison of shape descriptors for mice behavior recognition. In:Proceeding of Progress in Pattern Recognition, Image Analysis, Computer Vision, and Applications. Berlin Heidelberg, Germany:Springer, 2010. 370-377 [17] Wu Y, Lim J, Yang M H. Online object tracking:a benchmark. In:Proceedings of the 2013 IEEE Conference on Computer Vision and Pattern Recognition. Portland, OR, USA:IEEE, 2013. 2411-2418 [18] Smeulders A W M, Chu D M, Cucchiara R, Calderara S, Dehghan A, Shah M. Visual tracking:an experimental survey. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2013, 36(7):1442-1468 https://www.researchgate.net/profile/Rita_Cucchiara/publication/258924422_Visual_Tracking_An_Experimental_Survey/links/0f31752e12b97368b7000000.pdf [19] Li X, Hu W M, Shen C H, Zhang Z F, Dick A, Van Den Hengel A. A survey of appearance models in visual object tracking. ACM Transactions on Intelligent Systems and Technology, 2013, 4(4):Article No.58 https://arxiv.org/pdf/1303.4803.pdf [20] Thrun S, Burgard W, Fox D. Probabilistic Robotics. Cambridge, UK:MIT Press, 2006. [21] 李天成, 范红旗, 孙树栋.粒子滤波理论、方法及其在多目标跟踪中的应用.自动化学报, 2015, 41(12):1981-2002 http://www.aas.net.cn/CN/abstract/abstract18773.shtmlLi Tian-Cheng, Fan Hong-Qi, Sun Shu-Dong. Particle filtering:theory, approach, and application for multi-target tracking. Acta Automatica Sinica, 2015, 41(12):1981-2002 http://www.aas.net.cn/CN/abstract/abstract18773.shtml [22] Hartigan J A, Wong M A. Algorithm as 136:a K-means clustering algorithm. Journal of the Royal Statistical Society. Series C (Applied Statistics), 1979, 28(1):100-108 [23] Bruyninckx H. Open robot control software:the OROCOS project. In:Proceedings of the 2001 IEEE International Conference on Robotics and Automation. Seoul, South Korea:IEEE, 2001:2523-2528 期刊类型引用(1)
1. 张博，裴宇驰，黄慧敏. 基于多模板匹配的双模型自适应相关滤波跟踪算法. 太赫兹科学与电子信息学报. 2022(06): 618-625 . 
百度学术其他类型引用(6)
-
下载:
计量
- 文章访问数: 2048
- HTML全文浏览量: 441
- PDF下载量: 621
- 被引次数: 7
