CN101241601B - 一种图形处理的关节中心参数估计方法 - Google Patents

Abstract

本发明公开了一种图形处理的关节中心参数估计方法，其用于一通用计算机的运动捕获处理，包括以下步骤：对标志点进行刚性预处理；估计计算关节中心回归参数；从关节中心回归参数获得刚性点集；对刚性点集进行选择处理。本发明图形处理的关节中心参数估计方法由于采用了上述对标志点的刚性预处理，以及优化的关节中心初值估计方法，其在图形处理中消除了标志点的漂移，简化了处理过程，且提高了数据处理的准确性。

Description

一种图形处理的关节中心参数估计方法

技术领域

本发明涉及一种计算机的图形处理方法，尤其涉及的是一种模拟生物运动的显示处理的关节中心参数估计方法。

背景技术

运动捕获技术的发展历程可以追溯到20世纪70年代，当时，迪斯尼公司曾试图通过捕获演员的动作以改进动画制作效果。当计算机技术刚开始应用于动画制作时，纽约计算机图形技术实验室的Rebecca Allen就设计了一种光学装置，将演员的表演姿势投射在计算机屏幕上，作为动画制作的参考。从20世纪80年代开始，美国Biomechanics实验室、Simon Fraser大学、麻省理工学院等陆续开展了计算机人体运动捕获技术的研究。

此后，运动捕获技术吸引了越来越多研究人员和开发商的目光，并从试用性研究逐步走向了实用化。1988年，SGI公司开发了可捕获人头部运动和表情的运动捕获系统。目前，在发达国家，运动捕获已经进入了实用化阶段，已经有多家厂商相继推出了商品化的运动捕获设备，如Vicon、Polhemus、Sega Interactive、MAC、FilmBox、MotionAnalysis等。运动捕获技术的应用领域也远远超出了最开始的角色动画，并已经成功地应用于虚拟现实、游戏、人体工程学、运动模拟训练、生物力学研究等许多方面。

究竟什么是运动捕获？运动捕获技术领域的权威工程师AlbertoMenache在他的新书《Understanding Motion Capture for Computer Animationand Video Games》(对计算机动画和视频游戏的运动捕获的理解)中给出了一个准确的定义：“Motion capture is the process of recording a live motionevent and translating it into usable mathematical terms by tracking a number ofkey points in space over time and combining them to obtain a singlethree-dimensional representation of the performance.”(运动捕获就是记录一生物的活动事件并将之转换为可用的数字表示的过程，该过程通过追踪一系列的在一定时间和空间内的关键点的运动，并结合这些关键点以获取代表该行为的一单独三维表示)简单来说，运动捕获是将现场行为(LivePerformance)转换为数字化行为(Digital Performance)的技术。被捕获的主体可以是这个世界上任何具有行为能力的生物或物体；而记录下的关键点(Key Points)的位置应该在能最佳表示出主体不同肢体或部位运动的关键区域，如肢体中心位置或肢体相邻链接处。对于人体而言，这些关键点的位置就是人体的关节中心附近或人体骨骼突出的位置。在这些位置放置了一些传感器，比如反光球、红外光二级管或电磁传感器等能够往信号收集装置主动或被动传送信息的器件。这些器件在运动捕获领域中经常被称为标志点(Marker)，它们是运动捕获系统的重要组成部分之一。

然而，在现实生活当中，人们总是将角色动画(Character Animation)和运动捕获混为一谈，尽管两者有着本质的区别。运动捕获只是“收集”运动，而角色动画的核心是驱动角色运动的数据。为了获得一个运动的角色，我们需要进行运动捕获，并且把捕获到的运动数据“映射”到三维角色上去。然而，运动捕获记录下的只是那些可以描述运动的数据，通常是指标志点的三维空间坐标，而角色动画所需要的数据通常是指一组描述角色平动的平移数据和关节转动的角度数据。

从标志点的三维空间坐标数据到描述角色平动的平移数据和关节转动的角度数据的映射问题是一个很复杂的问题，也正是运动捕获数据处理技术研究的核心问题。它需要基于数字信号处理技术处理由于各种原因导致的三维空间坐标数据中的噪声；其次，它需要基于人体生理结构特性和运动特性，运用优化原理和计算机方法来生成人体骨骼运动。它的研究涉及到生理学、多刚体运动学、优化计算、图形学等多个领域，是一个多学科交叉的研究课题。

最近几年，计算机视觉领域对运动捕获数据处理问题进行了积极的探索，取得了很大的进展，也有了一些商业上的应用系统，比如Autodesk公司的MotionBuilder和Character Studio以及Vicon公司的Bodybuilder。这些软件没有技术文档可以参考，但是由于它们的数据处理结果中广泛存在膝关节和肘关节没伸展开和关节自由度被限制的情形，使得研究者们确信这些软件采用了逆运动学技术；另外这些软件有时会产生的头倒置的处理结果，使得研究者们猜测采用了一些诱导性的法则，例如：头的朝向被默认为始终朝上的。只有熟练的动画师才能使用这些软件，通过复杂的操作完成运动捕获数据的处理。因此，目前并没有完全令人满意的实用成果，尤其是离计算机动画、游戏和人体运动仿真等要求自动化准确恢复人体运动的目标还有相当大的距离。

另一方面，运动捕获技术在我国还是一个需要深入研究开发的领域，因为国内刚刚认识到它的作用，并且还只是刚刚参与进来。现有的两家厂商都是在机械复制国外的产品，在数据处理技术上遇到了瓶颈，无法提高数据的处理质量和简化运动捕获数据的处理操作流程。

由此可见，由于存在众多技术难点，已有的研究工作在精度、效率、稳定性等方面还难以满足实际应用中人们对高质量运动捕获数据的迫切需求。这就促使另辟蹊径，探讨运动捕获数据处理研究领域中的新方法，新思路。

现有技术的三维人体运动捕获技术是计算机视觉、计算机图形学以及虚拟现实等研究领域中一个备受关注的前沿方向，在计算机动画、人体运动仿真、生物运动力学分析、医疗康复等方面具有广阔的应用前景，不仅具有重要的研究意义，而且具有很好的应用价值。运动捕获技术的最大优点是能够捕获到物体真实运动的数据，由此生成的运动具有很强的真实感，并能合成更多复杂的运动。

准确的关节中心位置数据不仅可以增加虚拟人运动的真实感，而且还可以增强人体运动力学仿真分析的可信度，因此，如何利用标志点数据计算出高质量的关节中心位置数据是现有技术的研发主要问题——关节中心参数估计问题。

目前的关节中心估计方法从原理上来分主要有两大类：回归分析方法和函数分析方法，并且，基于函数分析方法的关节中心参数估计算法开始得到了高度的重视。然而在研究过程当中，研究者们更多的关注，针对不同情况发展一系列的不同目标函数，来进行关节中心参数的计算，而很少有人去关注对原始标志点数据的处理。

由于关节中心计算的目标函数中变量数目比较多(主要为关节中心的位置参数)，并且是非线性的，因此在做优化计算时，收敛性能很不好，主要表现在两个方面：(1)收敛速度比较慢；(2)容易在局部的极小值收敛。在优化计算当中，一个好的初值对于优化一个性质不好的目标函数而言，是至关重要的，它可以有效提高收敛速度，同时保证收敛到一个比较好的解。

在现有技术的运动捕获过程中，在表演者的一些特殊的肢体(如：头部，腰部)上，贴有一定数量(超过三个)的标志点，用来定位这些肢体的位姿，因此希望这些标志点在相应的肢体上构成一个刚体。然而，人体皮肤或衣服在表演者做出各种运动时，不可避免地会发生微小褶皱和位移，从而使得肢体上标志点的相对位置发生改变，最终会导致这些肢体位姿数据中含有大量的噪声，不仅会导致数据处理速度缓慢，而且会导致数据处理结果严重失真，甚至出现错误。如何有效的去除这些噪声，修正特定肢体上标志点的刚体性质就是标志点刚性修正问题。

估计两个对应点集的相对旋转、缩放和平移参数是一个重要的问题，并且被公认称为绝对定向问题。极小化式(1)这个最小二乘误差函数标志点刚性修正问题是比绝对定向问题更复杂的问题，其要解决的问题是：给定有点对应关系的多个点集，如何找到一个标准点集，并确定其他点集相对于它的朝向和位置，同时使得匹配误差最小。

$F(R,t,c)=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left|\right|y_i-c{\mathrm{Rx}}_i-t\left|\right|}^2---\left(8\right)$

解决绝对定向的迭代算法最早出现在二十世纪五十年代到六十年代之间，当时绝对定向技术主要应用在照相测量学(photogrammetry)当中，然后才是计算机视觉处理当中。与传统的非线性迭代相比，Thompson将绝对定向问题转换为线性方程组求解问题，可以说是绝对定向技术上的一大进步，尽管他当时只考虑了三个点构成一个点集的情况。在几个月以后，Schut依靠四元数极大地精简了相关线性方程组的推导过程。为了处理更多的点，Oswal和Balasubramanian发展了一种最小二乘方法，但是他们的方法只是找到一个最佳匹配的线性变换，因此需要做旋转矩阵正交化这样一个后处理。1973年，Sanso将旋转分量的求解归结为矩阵特征值分解，仍然使用了迭代方法，并不是解析(closed-form)解。1987年，Horn通过使用单位四元数(Unit Quaternion)代表旋转矩阵，改进了Sanso提出的方法。该方法既不需要解线性方程组，也不需要迭代更新，只需要解一个四次方程，是一种解析形式的解。由于Horn第一次给出绝对定向问题的解析解，因此他的方法引用最为广泛。

后经过十多年的发展，针对绝对定向问题，研究者们总共提出了四种不同的解析解法，这些方法的区别在于使用了不同的旋转变换表达形式和不同的目标函数。除了Horn的方法以外，1987年Arun提出了绝对定向问题的第二种解析方法，其中使用旋转矩阵表示旋转变换，并采用奇异值分解(Singular Value Decomposition)来进行求解。1988年，Horn又提出了第三种解析方法，其中正交矩阵被用来表示旋转变换，同样特征值分解被用来求解这个正交矩阵。在特殊的情况下，第二种方法和第三种方法会求得一个反射矩阵而非旋转矩阵，针对这种退化情形，Umeyama进行了修正，提出了一种比较完善的奇异值分解方法。Walker等人使用对偶四元数来表示旋转和平移变换，提出了第四种解析方法。

由于多视角数据配准(multiview registration)问题的核心是多个点集的配准，所以绝对定向技术被广泛应用到这个领域当中。其中，Williams和Bennamoun推广了Arun的方法，采用奇异值分解方法来计算最佳的配准参数。Krishnan等人基于牛顿方法提出了旋转群SO³上的一种迭代计算技术。

更一般的绝对定向问题是两个点集都受到噪声的污染，在这个假定下，首先就需要恢复出一个没有受到噪声污染的标准点集，然后再计算原始两个点集相对于标准点集的相似变换参数。针对这种情况，Goryn证明了：Arun提出的方法仍然可以达到一个较好的结果，然而并不能保证这个结果是一个最优解，甚至连局部最优解都保证不了。Ramos等人则提出了一种混合最小二乘方法，但是他们的方法得到的结果并不能保证是一个旋转矩阵，需要做旋转矩阵的正交化这样一个后处理。因此，上述两项研究工作并未能完善解决这种无标准点集的绝对定向问题。

实际当中碰到的问题经常是，点集的个数远远大于二，在这种情况下如何求取点集之间的相似变换参数，即对于多点集最小二乘拟合问题，在运动捕获数据处理，皮肤变形估计，图像注册(Image Registration)以及模式识别当中有广泛的应用。wen等人给出了基于绝对定向技术和梯度下降法的多点集最小二乘问题的迭代解法。

关节中心估计问题一直是运动捕获技术和生物力学分析技术等研究领域中的一个热点，至今仍然有不少国内外学者在该问题上开展研究工作。目前，关节中心的估计方法主要分为两大类：(1)基于回归分析的关节中心估计方法；(2)基于函数分析的关节中心估计方法。

基于回归参数的方法很容易实现，速度也非常快，实用于实时人体运动捕获，是目前应用最广的一种关节中心估计算法。但是回归分析和测量设备容易引入误差，所以这种标准化方法的精度是有限的，误差可以达到25-30毫米。

基于函数分析的关节中心估计方法的理论基础是：相邻两个刚体的公共点就是关节中心。基于这个理论，O′Brien等人针对电磁学运动捕获数据，采用奇异值分解的方法估计了人体关节中心。在这个方法中，由于电磁学运动捕获数据可以提供每个肢体的局部坐标系，因此可以很方便的实施该求解方法。基于运动学约束，Schwartz and Rozumalski提出了一个改进的方法，能够更加合理的计算关节中心和旋转轴。

如上所述，基于函数分析的关节中心计算方法已经积累不少的经验，但是由于关节中心的参数化坐标系方式不一样，并且目标函数的构造不一样，根据关节中心的参数化坐标和目标函数的构造，可以对基于函数分析的关节中心计算方法进行更细致的分类。

在现有技术的种种图形处理方法和运动捕获技术中，在演员的一些特殊的肢体(头部，腰部)上贴有一定数量(超过三个)的标志点，一般希望这些肢体上的标志点构成一个刚体，从而用来定位这些肢体的位姿。然而，人体皮肤或衣服在演员做出各种运动时，不可避免地会发生微小褶皱和位移，从而使得肢体上标志点的相对位置会发生改变，最终导致这些肢体位姿数据含有大量噪声。如何有效的去除这些噪声，修正特定肢体上标志点的刚体性质就是本发明要解决的第一个主要问题——标志点数据刚性修正问题。

准确的关节中心位置数据不仅可以增加虚拟人运动的真实感，而且还可以增强人体运动力学仿真分析的可信度，因此，如何利用标志点数据计算出高质量的关节中心数据就是本发明要解决的第二个主要问题——关节中心参数估计问题。

发明内容

本发明的目的在于提供一种图形处理的关节中心参数估计方法，针对现有技术的缺陷，采用最小二乘匹配对于标志点的刚性预处理，以及快速基于旋转几何和非线性优化的关节中心初值估计方法，实现准确的骨骼之关节中心参数估计和处理。

本发明的技术方案如下：

一种图形处理的关节中心参数估计方法，其用于一通用计算机的运动捕获处理，包括以下步骤：

A、对标志点进行刚性预处理；所述的步骤A包括：

A1、输入点集{X^j}_j＝1 ^k；

A2、计算无噪声的点集

相似变换参数{R_j}_j＝1 ^k，{t_j}_j＝1 ^k和{c_j}_j＝1 ^k；

A3、i从1到n，同时j从1到k获得输出修正点集{Y^j}_j＝1 ^k的元素 $y_i^j=c_j{R}_j{z}_i+t_j;$

B、估计计算关节中心回归参数；所述的步骤B包括：

B1、输入点集{P_i ^j}_i＝1 ³和{C_i ^j}_i＝1 ³；

B2、采用非线性优化极小化误差函数 $\mathrm{ERROR}=\underset{i=1}{\overset{F}{\mathrm{\Σ}}}{({\left|\right|J_i-C_1^i\left|\right|}^2-r^2)}^2;$

其中，F为采集数据的帧数，ⁱJ为关节中心，ⁱC₁为子刚体上的标志点，r为回归坐标参数；

B3、输出回归参数{a，b，c}；

C、从关节中心回归参数获得刚性点集；所述的步骤C包括：

C1、更新点集{P_i ^j}_i＝1 ³和{C_i ^j}_i＝1 ³，去除噪声；

C2、计算矩阵A，B：采用相邻的父子刚体上各三个点，ⁱP₁、ⁱP₂、ⁱP₃是父刚体上的标志点，角标i表示帧数，ⁱC₁、ⁱC₂、ⁱC₃是子刚体上的标志点，关节中心ⁱJ相对于相邻刚体上的标志点始终保持不变，满足如下关系式：

ⁱJ＝[ⁱP₁-ⁱP₂ ⁱP₁-ⁱP₃ (ⁱP₁-ⁱP₂)×(ⁱP₁-ⁱP₃) ⁱP₁]·[a b c 1]^t    (9)

ⁱJ＝[ⁱC₁-ⁱC₂ ⁱC₁-ⁱC₃ (ⁱC₁-ⁱC₂)×(ⁱC₁-ⁱC₃) ⁱC₁]·[d e f 1]^t    (10)

其中，{a，b，c，d，e，f}为回归坐标参数；

两式联立可得到：

[ⁱP₁-ⁱP₂ ⁱP₁-ⁱP₃ (ⁱP₁-ⁱP₂)×(ⁱP₁-ⁱP₃)  ⁱC₁-ⁱC₂ ⁱC₁-ⁱC₂ (ⁱC₁-ⁱC₂)×(ⁱC₁-ⁱC₃)]·

[a  b  c  -d  -e  -f]^t＝-(ⁱP₁-ⁱC₁)                          (11)

对于连续采集的F帧数据，得到如下方程组：

AX＝B                                                       (12)

$A=\left[\begin{array}{cccccc}P_1^1-P_2^1& P_1^1-P_3^1& (P_1^1-P_2^1)\×(P_1^1-P_3^1)& C_1^1-C_2^1& C_1^1-C_2^1& (C_1^1-C_2^1)\×(C_1^1-C_3^1)\\ P_1^2-P_2^2& P_1^2-P_3^2& (P_1^2-P_2^2)\×(P_1^2-P_3^2)& C_1^2-C_2^2& C_1^2-C_2^2& (C_1^2-C_2^2)\×(C_1^2-C_3^2)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ P_1^F-P_2^F& P_1^F-P_3^F& (P_1^F-P_2^F)\×(P_1^F-P_3^F)& C_1^F-C_2^F& C_1^F-C_2^F& (C_1^F-C_2^F)\×(C_1^F-C_3^F)\end{array}\right]---\left(13\right)$

X＝[a b c -d -e -f]^t

B＝[-(¹P₁-¹C₁)-(²P₁-²C₁)…-(^FP₁-^FC₁)]^t                      (14)

C3、计算矩阵X：当采集数据的帧数F超过3时，方程组变成超定方程组，可以最小二乘解为：

X＝(A^tA)^-1A^tB                                               (15)

通过式(14)或(15)计算出每一帧关节中心的位置；

C4、输出回归参数{a，b，c，d，e，f}；

D、对刚性点集进行选择处理；所述的步骤D包括：

D1、输入点集

处理点集

得到{z_i}_i＝1 ^K，变换参数{R_j}_i＝1 ^F、{t_i}_i＝1 ^F，{cⁱ}_i＝1 ^F；

D2、i从1到K计算 $\{f_i=\frac{1}{F}\underset{j=1}{\overset{F}{\mathrm{\Σ}}}{\left|\right|y_i^j-c_j{R}_j{z}_i-t_j\left|\right|}^2\};$

D3、按递减顺序给f_i排序成f_i0≤f_i1≤...≤f_ik-1；

D4、按序号i1，i2，...，ik-1，从中选择与z_i0，z_i1不共线的第三点z_it，输出i0，i1，it。

本发明所提供的一种图形处理的关节中心参数估计方法，由于采用了上述对标志点的刚性预处理，以及优化的关节中心初值估计方法，其在图形处理中消除了标志点的漂移，简化了处理过程，且提高了数据处理的准确性。

附图说明

图1为本发明方法的处理过程一的示例图；

图2为本发明方法的处理过程二的示例图；

图3为本发明方法的标志点在手臂上的设置示意图；

图4(a)和图4(b)为本发明方法的刚性点集选择算法统计结果示意图；

图5(a1)-(a3)为大范围连续数据下关节中心算法性能在不同比较噪声比下的对比示意图；

图5(b1)-(b3)为大范围非连续数据下关节中心算法性能比较在不同噪声比下的对比示意图；

图6为本发明方法的标志点在人体上的分布示意图；

图7(a)-(d)为对人体的走、跳、ROM运动以及跑的关节中心计算处理的示意图；

图8为本发明方法与现有技术对人体某个肢体的运动计算处理结果的效果比对示意图；

图9为本发明方法较佳实施例与现有技术对四肢长度数据的标准方差对比效果示意图；

图10为本发明方法较佳实施例与现有技术对四肢长度数据的均值对比示意图。

具体实施方式

以下对本发明的较佳实施例加以详细说明。

本发明的关节中心参数估计方法中，由于人体皮肤或者衣服的微小褶皱和位移，而引起肢体上标志点的相对位置的改变，而导致这些肢体位姿数据中含有大量噪声。为有效去除各种噪声，提高标志点的数据质量，本发明方法采用了基于成熟的绝对定向技术和采用梯度下降逐步求解的标志点刚性预处理算法，极小化下面这个最小二乘误差函数：

$F(Z,R,t,c)=\frac{1}{k}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left|\right|y_i^j-c_j{R}_j{z}_i-t_j\left|\right|}^2---\left(16\right)$

其具体的处理过程如下：

一、标志点的刚性预处理过程如下：

在本发明的图形处理计算机系统中输入：点集{X^j}_j＝1 ^k；输出：修正点集{Y^j}_j＝1 ^k；

其步骤包括：(A1)计算无噪声的点集

相似变换参数{R_j}_j＝1 ^k，{t_j}_j＝1 ^k和{c_j}_j＝1 ^k；(A2)For i＝1 to n

$\mathrm{For\; j}=1\mathrm{to\; k}\{y_i^j=c_j{R}_j{z}_i+t_j;\}.$

由于本发明方法的关节中心计算目标函数中变量数目比较多，主要为关节中心的位置参数，并且是非线性的，因此在做优化计算时，收敛性能很不好，其主要表现在两个方面：(1)收敛速度比较慢，意味着运算速度太慢；(2)容易在局部的极小值收敛，意味着会得到错误的关节中心估计参数值。在本发明方法优化处理过程当中，一个好的初值对于优化一个性质不好的目标函数而言，是至关重要的，它可以有效提高收敛速度，同时保证收敛到一个比较好的解。因此，本发明图形处理的关节中心参数估计方法即可获得较好的初值。

在人体运动的人体结构上，相邻的骨骼可以看成通过关节连接的刚体，在图形处理中，第i帧相邻两个刚体中，选其中为一父刚体，另一相邻刚体为子刚体，则父、子刚体中一个需要有三个点，另外一个则需要少于三个点。由于对于旋转中心而言，父、子刚体是对称的，所以本发明方法假定父刚体有三个点，而子刚体只有一个点，如图1所示的。其中ⁱP₁、ⁱP₂、ⁱP₃是父刚体上的标志点，角标i表示帧数，ⁱC₁是子刚体上的标志点，由旋转几何的知识可知，关节中心ⁱJ相对于相邻刚体上的标志点始终保持不变，由父刚体上的三点有如下关系式：

ⁱJ＝[ⁱP₁-ⁱP₂ ⁱP₁-ⁱP₃ (ⁱP₁-ⁱP₂)×(ⁱP₁-ⁱP₃) ⁱP₁]·[a b c 1]^t    (17)

其中，a、b、c为回归坐标参数。同时，由于标志点与关节中心的距离保持不变，由子刚体上的一点有：

||ⁱJ-ⁱC₁||²＝r²                                             (18)

两式联立可得到一个误差函数：

Error_i＝(||ⁱJ-ⁱC₁||²-r²)²                                   (19)

对于连续采集的F帧数据，可以得到如下总的误差函数：

$\mathrm{ERROR}=\underset{i=1}{\overset{F}{\mathrm{\Σ}}}{\mathrm{Error}}_i=\underset{i=1}{\overset{F}{\mathrm{\Σ}}}{({\left|\right|J_i-C_1^i\left|\right|}^2-r^2)}^2---\left(20\right)$

采用非线性优化对上式极小化，就可以获取参数回归坐标参数a、b、c、r。然后通过通用计算装置，例如计算机，计算出每一帧关节中心的位置。如果子刚体上有两个点，只需要在式(10)中增加类似的一个误差项就可以了。

以下说明本发明方法中采用关节中心回归参数估计的处理过程：

输入：点集{P_i ^j}_i＝1 ³和{C_i ^j}_i＝1 ³。

输出：回归参数{a，b，c}

其处理步骤包括：(B1)根据(A1)和(A2)的步骤更新点集{P_i ^j}_i＝1 ³；

(B2)采用非线性优化极小化误差函数 $\mathrm{ERROR}=\underset{i=1}{\overset{F}{\mathrm{\Σ}}}{({\left|\right|J_i-C_1^i\left|\right|}^2-r^2)}^2;$

(B3)输出回归参数{a，b，c}。

此时采用相邻的父子刚体上各三个点，如图2所示：其中ⁱP₁、ⁱP₂、ⁱP₃是父刚体上的标志点，角标i表示帧数，ⁱC₁、ⁱC₂、ⁱC₃是子刚体上的标志点，由旋转几何的知识可知，关节中心ⁱJ相对于相邻刚体上的标志点始终保持不变，满足如下关系式：

ⁱJ＝[ⁱP₁-ⁱP₂ ⁱP₁-ⁱP₃ (ⁱP₁-ⁱP₂)×(ⁱP₁-ⁱP₃) ⁱP₁]·[a b c 1]^t    (21)

ⁱJ＝[ⁱC₁-ⁱC₂ ⁱC₁-ⁱC₃ (ⁱC₁-ⁱC₂)×(ⁱC₁-ⁱC₃) ⁱC₁]·[d e f 1]^t    (22)

其中，{a，b，c，d，e，f}为回归坐标参数。

两式联立可得到：

[ⁱP₁-ⁱP₂ ⁱP₁-ⁱP₃ (ⁱP₁-ⁱP₂)×(ⁱP₁-ⁱP₃) ⁱC₁-ⁱC₂ ⁱC₁-ⁱC₂ (ⁱC₁-ⁱC₂)×(ⁱC₁-ⁱC₃)]·

[a b c -d -e -f]^t＝-(ⁱP₁-ⁱC₁)                               (23)

对于连续采集的F帧数据，本发明方法可以得到如下方程组：

AX＝B                                                       (24)

其中：

$A=\left[\begin{array}{cccccc}P_1^1-P_2^1& P_1^1-P_3^1& (P_1^1-P_2^1)\×(P_1^1-P_3^1)& C_1^1-C_2^1& C_1^1-C_2^1& (C_1^1-C_2^1)\×(C_1^1-C_3^1)\\ P_1^2-P_2^2& P_1^2-P_3^2& (P_1^2-P_2^2)\×(P_1^2-P_3^2)& C_1^2-C_2^2& C_1^2-C_2^2& (C_1^2-C_2^2)\×(C_1^2-C_3^2)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ P_1^F-P_2^F& P_1^F-P_3^F& (P_1^F-P_2^F)\×(P_1^F-P_3^F)& C_1^F-C_2^F& C_1^F-C_2^F& (C_1^F-C_2^F)\×(C_1^F-C_3^F)\end{array}\right]---\left(25\right)$

X＝[a b c -d -e -f]^t

B＝[-(¹P₁-¹C₁)-(²P₁-²C₁)…-(^FP₁-^FC₁)]^t                      (26)

当采集数据的帧数F超过3时，方程组变成超定方程组，可以最小二乘解为：

X＝(A^tA)^-1A^tB                                               (27)

通过式(16)或(17)可以计算出每一帧关节中心的位置。

本发明方法中的关节中心回归参数估计得到刚性点集的处理过程包括：

输入：点集{P_i ^j}_i＝1 ³和{C_i ^j}_i＝1 ³。

输出：回归参数{a，b，c，d，e，f}

步骤：(C1)更新点集{P_i ^j}_i＝1 ³和{C_i ^j}_i＝1 ³，去除噪声。

(C2)计算矩阵A，B。

(C3)计算矩阵X；

(C4)输出回归参数{a，b，c，d，e，f}。

然后本发明方法需要对计算出的点集进行选择：

当刚体上的点超过3个时，假定第i帧父刚体为{ⁱP_j}_j＝1 ^N(N＞3)，。实际上，在刚性条件下，每个刚体上任意三个不共线的标志点对关节中心的位置起到的作用是等效的，因此，只需要在刚体上选取三个不共线的标志点就可以了。本发明方法期望能够选择出最稳定的三个点，因此，参考前述刚性算法处理过程，本发明选择修正前后误差最小的不共线的三个标志点。

本发明方法中的刚性点集选择处理过程包括：

输入：点集

输出：选择点的标号t₁，t₂，t₃。

具体的处理步骤包括：(D1)参考步骤(A1)(A2)的处理，处理点集

得到{z_i}_i＝1 ^K，变换参数{Rⁱ}_i＝1 ^F、{tⁱ}_i＝1 ^F，{cⁱ}_i＝1 ^F；

(D2) $\mathrm{for\; i}=1\mathrm{to\; K\; do}\{f_i=\frac{1}{F}\underset{j=1}{\overset{F}{\mathrm{\Σ}}}{\left|\right|y_i^j-c_j{R}_j{z}_i-t_j\left|\right|}^2\};$

(D3)按递减顺序给f_i排序成f_i0≤f_i1≤...≤f_ik-1；

(D4)按序号i1，i2，...，ik-1，从中选择与z_i0，z_i1不共线的第三点z_it，输出i0，i1，it.

当父、子刚体的三个点都选出来以后，就可以用步骤(C1)到(C4)的处理进行进一步的关节中心求解计算。

本发明方法可以在通用计算机使用matlab 7.0实现了上述各处理过程，与已有的一些比较成熟的关节中心参数估计算法Silaghi1998，OBrein00，UdLa02，KOF05，进行了对比测试，并且在一台个人电脑(Pentium IV，2.8GHz，内存1.0GB)上进行了性能测试。

如图3所示的，本发明采用手臂做为部分示例说明，比较本发明的标志点刚性预处理方式、绝对定向算法(OA)、基于小波滤波器的方法(Wden)和线性滤波器方法(LSI)。其中Wden方法直接使用matlab函数wden实现，参数为(x，′heursure′，′s′，′one′，3，′sym8′)，即针对基本的高斯噪声采用8阶的Daubechies小波对信号进行三次小波分解后，使用混合软阈值进行噪声去除。线性滤波器方法(LSI)采用5阶的滤波器(1/16，4/16，6/16，4/16，1/16)进行卷积滤波过滤处理。

为了评价过滤之后和对关节中心计算质量的改进程度，定义相对误差下降比例EDR：

做两个实验：第一个实验中，标志点数据有大范围连续和小范围连续两组。第二个实验中，标志点数据由随机生成的大范围离散和小范围离散两组数据构成。在不同级别的噪声下，每一组实验进行100次试验，记录下平均响应作为该组实验的结果。

首先通过沿着特定曲线取旋转参数值使子刚体上的点绕关节中心J旋转，然后通过沿着特定的曲线同时旋转和平移父刚体和子刚体上的点，再加上服从高斯分布的随机噪声之后形成大范围离散数据的测试点集。通过控制其中旋转参数的取值范围：[-5°，5°]×[-5°，5°]×[-5°，5°]，生成小范围连续测试数据。

随机选择的旋转参数使子刚体上的点绕关节中心J旋转，然后通过随机选择的旋转参数和平移参数同时旋转和平移父刚体和子刚体上的点，再加上服从高斯分布的随机噪声之后形成测试点集其中参数旋转和平移参数的取值范围有两组分别为：I.[-5°，5°]×[-5°，5°]×[-5°，5°]和[-1，1]和II.[-180°，180°]×[-90°，90°]×[-90°，90°]，[-100，100]。

实验结果发现，在数据噪声水平为5％的前提下，Si1aghi98、OBrein00、UdLa02、KOF05四种关节中心算法的相对误差下降比例，在大范围连续测试数据集下分别可以达到28.62％、7.44％、2.30％、17.33％；在小范围连续测试数据集下分别可以达到14.57％、13.42％、13.57％、14.23％；在大范围离散数据集下分别可以达到35.53％、6.48％、3.01％、35.08％；在小范围离散数据集下分别可以达到15.57％、11.68％、11.70％、13.26％；总的平均相对误差下降比例分别可以达到23.57％、9.76％、7.64％、19.97％，由此可见基于标志点数据刚性修正算法的数据预处理步骤可以显著提高关节中心参数估计算法的质量。

参照图3所示的父子刚体模型，从每个刚体上的五个点中选出三个性质好的点来参与关节中心计算。对于标志点受到同等水平噪声污染和不同水平噪声做试验测试，使用步骤(D1)到(D4)实现了刚性点集选择算法，并使用步骤(C1)到(C4)实现了对选择出来的点集进行关节中心的计算。

在第一组实验中，每一次实验中，每一个标志点加的噪声都是同样水平的，从五个点中选三个点一共是

种选法，父子刚体选择标志点的方案有

种。在相同的噪声水平下，比较本发明方法计算关节中心和原始仿真关节中心平均距离最优性统计结果和相应的方差统计结果，发现没有一个绝对占优的父子刚体标志点选择组合，每种组合计算的结果非常接近，因此任意选择一组父子刚体标志点选择组合进行关节中心计算都是可行的。

造成这种现象的原因是事先模拟的数据点中，每一个标志点上加的噪声都是同等水平的，如果每个标志点上加入的噪声不一样，也就不能任意选择一组父子刚体标志点选择组合进行关节中心计算。如果在标志点上加入的噪声是同等水平的情况下，都能选择一组较好的组合，那么可以预期本发明方法在不同噪声水平情况下的效果会更好。在图4中给出了对应的统计结果，图4(a)示出的是参数II下计算关节中心和原始仿真关节中心平均距离排名图，图4(b)示出的是相应的方差统计结果，显然步骤(D1)到(D4)可以提出一个较好的父子刚体标志点组合选择方案，需要注意的是这种方案选择的组合是变化的。

从上述讨论，可以知道：对于第二类方法，搜索目标函数的最小值时很容易就陷入到局部极值点，因此如果想快速收敛到一个比较好的解，就必须有一个好的初值，步骤(C1)到(C4)可以快速提供这样一个初值。本发明方法实现了基于目标函数{f_i}_i＝1 ⁶和f₃ ^*(α＝0.01)的关节求解算法，经过比较，发现基于目标函数f₄的关节中心算法性能是最好的。

具体实验如下进行，使用四份测试数据，分别进行关节中心求解计算，统计结果如下。其中，关节中心误差评价函数F如下定义：

其中，T为试验次数，^tm_k为优化计算所得的第t次试验第k帧关节中心，^tm_k ^*为第t次试验第k帧原始仿真关节中心。

如图5(a1)-(a3)以及图5(b1)-(b3)所示，大范围连续数据和大范围非连续数据下的实验结果表明：基于目标函数f₄的关节中心算法恢复的关节中心比其他的方法计算所得的关节中心的精度要高，稳定性要好(方差最小)，这一点随着噪声水平的增加而更加明显。基于目标函数f₄的优化计算恢复的关节中心要比基于目标函数f₃ ^*(α＝0.01)的优化计算恢复的关节中心精度要高、稳定性要好，这一点随着噪声水平的增加也更加明显。

下面就运动捕获技术数据的处理说明对本发明技术的应用：

本发明方法使用MATLAB实现使用稀疏矩阵的关节中心位置参数化方法来解决链状骨骼匹配的方法，并使用OpenGL图形库建立一个三维虚拟场景。使用VICON运动捕获设备，捕获了几类人体运动，包括走(76帧)、跑(47帧)和跳(105帧)和一个特殊的关节活动运动(Range OfMotion)(1085帧)。在这个关节活动运动中，演员被要求尽可能的伸展肢体和活动关节。运动捕获使用了12个MX-40型号的摄像机，采样频率是30Hz，演员的身上贴有59个标志点，保证四肢的肢体上都有三个标志点，如图6所示。

本发明方法分别计算了人体四肢12个主要的关节中心，如图7(a)-图7(d)所示的，显示了处理ROM运动、走、跳、跑运动捕获数据的处理结果示意图。

同时用f₁、f₂、f₃ ^*(α＝0.01)、f₄所对应的关节中心算法处理了上述数据，然后统计了相应的骨骼长度，图8显示了左大腿长度的计算结果。从方差的角度来考察，该图说明，第二类方法的关节中心(f₃ ^*(α＝0.01)、f₄)计算恢复的关节中心比第一类方法(f₁、f₂)计算所得的关节中心的稳定性要好。从均值的角度来考察，该图说明，目标函数f₃ ^*(α＝0.01)中额外添加的约束项，会影响到关节中心计算的精度，而基于目标函数f₄的优化计算恢复的关节中心要比基于目标函数f₃ ^*(α＝0.01)的优化计算恢复的关节中心精度要高。

在图9和图10中，对比四肢的骨骼长度数据的标准方差和均值，结果表明：与其它三种关节中心估计算法相比，基于目标函数f₄的关节中心求解算法所对应的标准方差是最小的，因此它是一种结果高稳定性的方法。与基于目标函数f₃ ^*(α＝0.01)的关节中心估计算法相比，基于目标函数f₄的关节中心求解算法所对应的关节中心数据与其它两种方法的结果更为接近，因此它是一种高精度的方法。

通过仿真数值试验验证了基于标志点数据刚性修正算法的数据预处理步骤可以显著提高关节中心参数估计算法的质量；在此基础上，本发明方法提出了一种高精度、高稳定性的基于函数分析的关节中心算法，这种算法优点在于提供了一种解析的关节中心初值计算方法，并且由这种方法计算得到的初值参与特定的目标函数优化计算后所得的关节中心，与以往的关节中心算法计算所得到的关节中心相比，具有更高的精度、更高的稳定性。通过大量的合成数据以及实际运动捕获数据处理实验，验证了此算法的有效性和实用性。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (1)

1.一种图形处理的关节中心参数估计方法，其用于一通用计算机的运动捕获处理，包括以下步骤：

A、对标志点进行刚性预处理；所述的步骤A包括：

A1、输入点集{X^j}_j＝1 ^k；

A2、计算无噪声的点集

相似变换参数{R_j}_j＝1 ^k，{t_j}_j＝1 ^k和{c_j}_j＝1 ^k；

A3、i从1到n，同时j从1到k获得输出修正点集{Y^j}_j＝1 ^k的元素

$y_i^j=c_j{R}_j{z}_i+t_j;$

B、估计计算关节中心回归参数；所述的步骤B包括：

B1、输入点集{P_i ^j}_i＝1 ³和{C_i ^j}_i＝1 ³；

B2、采用非线性优化极小化误差函数 $\mathrm{ERROR}=\underset{i=1}{\overset{F}{\mathrm{\Σ}}}{({\left|\right|J_i-C_1^i\left|\right|}^2-r^2)}^2;$

其中，F为采集数据的帧数，ⁱJ为关节中心，ⁱC₁为子刚体上的标志点，r为回归坐标参数；

B3、输出回归参数{a，b，c}；

C、从关节中心回归参数获得刚性点集；所述的步骤C包括：

C1、更新点集{P_i ^j}_i＝1 ³和{C_i ^j}_i＝1 ³，去除噪声；

C2、计算矩阵A，B：采用相邻的父子刚体上各三个点，ⁱP₁、ⁱP₂、ⁱP₃是父刚体上的标志点，角标i表示帧数，ⁱC₁、ⁱC₂、ⁱC₃是子刚体上的标志点，关节中心ⁱJ相对于相邻刚体上的标志点始终保持不变，满足如下关系式：

ⁱJ＝[ⁱP₁-ⁱP₂ ⁱP₁-ⁱP₃(ⁱP₁-ⁱP₂)×(ⁱP₁-ⁱP₃)ⁱP₁]·[a b c 1]^t      (1)

ⁱJ＝[ⁱC₁-ⁱC₂ ⁱC₁-ⁱC₃(ⁱC₁-ⁱC₂)×(ⁱC₁-ⁱC₃)ⁱC₁]·[d e f 1]^t      (2)

其中，{a，b，c，d，e，f}为回归坐标参数；

两式联立可得到：

[ⁱP₁-ⁱP₂ ⁱP₁-ⁱP₃(ⁱP₁-ⁱP₂)×(ⁱP₁-ⁱP₃)ⁱC₁-ⁱC₂ ⁱC₁-ⁱC₂(ⁱC₁-ⁱC₂)×(ⁱC₁-ⁱC₃)]·[a b c-d-e-f]^t＝-(ⁱP₁-ⁱC₁)                                (3)

对于连续采集的F帧数据，得到如下方程组：

AX＝B                                                 (4)

其中：

$A=\left[\begin{array}{cccccc}P_1^1-P_2^1& P_1^1-P_3^1& (P_1^1-P_2^1)\×(P_1^1-P_3^1)& C_1^1-C_2^1& C_1^1-C_2^1& (C_1^1-C_2^1)\×(C_1^1-C_3^1)\\ P_1^2-P_2^2& P_1^2-P_3^2& (P_1^2-P_2^2)\×(P_1^2-P_3^2)& C_1^2-C_2^2& C_1^2-C_2^2& (C_1^2-C_2^2)\×(C_1^2-C_3^2)\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ P_1^F-P_2^F& P_1^F-P_3^F& (P_1^F-P_2^F)\×(P_1^F-P_3^F)& C_1^F-C_2^F& C_1^F-C_2^F& (C_1^F-C_2^F)\×(C_1^F-C_3^F)\end{array}\right]---\left(5\right)$

X＝[a b c-d-e-f]^t

B＝[-(¹P₁-¹C₁)-(²P₁-²C₁)…-(^FP₁-^FC₁)]^t                (6)

C3、计算矩阵X：当采集数据的帧数F超过3时，方程组变成超定方程组，可以最小二乘解为：

X＝(A^tA)^-1A^tB                                         (7)

通过式(6)或(7)计算出每一帧关节中心的位置；

C4、输出回归参数{a，b，c，d，e，f}；

D、对刚性点集进行选择处理；所述的步骤D包括：

D1、输入点集

，处理点集

得到{z_i}_i＝1 ^K，变换参数{Rⁱ}_i＝1 ^F、{tⁱ}_i＝1 ^F，{cⁱ}_i＝1 ^F；

D2、i从1到K计算 $\{f_i=\frac{1}{F}\underset{j=1}{\overset{F}{\mathrm{\Σ}}}{\left|\right|y_i^j-c_j{R}_j{z}_i-t_j\left|\right|}^2\};$

D3、按递减顺序给f_i排序成f_i0≤f_i1≤...≤f_ik-1；

D4、按序号i1，i2，...，ik-1，从中选择与z_i0，z_i1不共线的第三点z_it，输出i0，i1，it。

Images