引擎设计跟踪(九.14.2 final) Inverse Kinematics: CCD 在Blade中的实现

因为工作忙, 好久没有记笔记了, 但是有时候发现还得翻以前的笔记去看, 所以还是尽量记下来备忘.

 关于IK, 读了一些paper, 觉得之前翻译的那篇, welman的paper (http://graphics.ucsd.edu/courses/cse169_w04/welman.pdf  摘译:http://www.cnblogs.com/crazii/p/4662199.html) 非常有用, 入门必读. 入门了以后就可以结合工程来拓展了.

先贴一下CCD里面一个关节的分析:

当Pic的方向和Pid重合时, 末端器离目标的距离最近, 所以把Pic绕着旋转轴旋转Φ度就可以. 当然旋转轴有方向,Φ也有方向(顺时针/逆时针).

如果取转轴axis  = Pic x Pid, 如果Pic和Pid已经单位化, 那么旋转角度等于asin(|axis|).

实际中可能会先判断两个向量是否已经同向, 如果方向不一致再旋转. 代码可以简单表示如下:

Pic.normalize();
Pid.normalize();
scalar cosAngle = Pic.dotProduct(Pid);
if( cosAngle < 1.0f ) //whether in same direction
{
    Vector3 axis = Pic.crossProduct(Pid);
    scalar angle = std::acos(cosAngle);
    Quaternion rotation(axis, angle);
    //rotate the joint using rotation
}

这里跟welman paper里的思路相同, 但是有细微的不同:

• 这里直接使用了3D的vector math来旋转向量, 而welman用了几何分析来建模, 之后就用代数简约并用导数求极值.
• 这里考虑到关节的多个DOF(Dimensions Of Freedom), axis是根据Pic和Pid两个方向叉积直接求得, 故axis可以是任意方向, 而原文中使用的是已知的固定转轴.
• 原文中的CCD没有奇异性, 因为用的固定转轴, 但是这里的方法是有奇异方向的: 当Pic和Pid一开始就方向相同或相反方向的时候, (跟原文中雅克比转置的奇异方向一样), 已经不需要旋转. 况且Pic x Pid是0向量, 转轴丢失, 产生了Gimbal lock. 这个问题后面再分析.

 上面是对CCD问题的基本算法, 下面记录实际使用时的一些问题:

约束(constraints)

考虑到DOF和旋转角度的限制, 需要给每个关节的旋转角度加上最大值和最小值. 这个些限制是关节点的局部限制, 并且与CCD旋转时的那个旋转过程量的转轴无关, 可以使用用局部旋转的欧拉角.

这就需要上面CCD迭代旋转关节时, 最好在关节的局部空间进行. 当然我也试过在模型空间(这里可以等同于世界空间)计算CCD, 不过应用约束的时候必须转换到局部空间.

约束角度的建立, 可以根据现实中关节的角度来定义, 比如"胳膊肘不能往外拐", 来限制关节的旋转. 需要注意的是, constraint的角度定义要跟artist建模时的空间一致.
比如建模时, 模型正面朝向z+, 那么膝盖的转轴就是x, 对应的欧拉角分量为pitch. 当模型正面朝向x+, 那么膝盖的转轴就变成z, 对应roll. 所以如果要配置约束角, 要跟实际美术的规范一致.

不过Blade角度约束是这样得来的: 分析原始FK动画, 得到所有关节的活动范围, 把它作为约束.
这是在快速阅读某个paper时发现的方法, 不过因为看的paper太多, 这个方法也只是在文中一句带过, 所以一不小心就可能没注意到. 还有一个Inverse Inverse Kinematics的方法, 也是分析FK来求解IK的, 不过没有读.
这个方法非常的简单, 缺点是需要有原始FK的复杂完整动画才有效. 比如一个简单站立动画, 腿部从来没有弯曲过, 膝盖没有旋转, 那么生成的约束范围就太小, 导致IK不可能出现弯腿.

对于原始FK动画的分析, blade是把它放在动画导入/导出时做的, 也就是在生成动画文件的时候, 顺便提取了所有关键中的数据, 并保存在骨骼动画文件中.
blade中constraints的定义如下:

    typedef struct IKConstraints
    {
        fp32    mMinX;
        fp32    mMaxX;
        fp32    mMinY;
        fp32    mMaxY;
        fp32    mMinZ;
        fp32    mMaxZ;
        inline IKConstraints()
        {
            mMinX = mMinY = mMinZ = FLT_MAX;
            mMaxX = mMaxY = mMaxZ = -FLT_MAX;
        }
        inline IKConstraints(fp32 minx, fp32 maxx, fp32 miny, fp32 maxy, fp32 minz, fp32 maxz)
        {
            mMinX = minx; mMaxX = maxx;
            mMinY = miny; mMaxY = maxy;
            mMinZ = minz; mMaxZ = maxz;
        }

        inline void merge(const SIKConstraints& rhs)
        {
            mMinX = std::min(mMinX, rhs.mMinX);
            mMaxX = std::max(mMaxX, rhs.mMaxX);
            mMinY = std::min(mMinY, rhs.mMinY);
            mMaxY = std::max(mMaxY, rhs.mMaxY);
            mMinZ = std::min(mMinZ, rhs.mMinZ);
            mMaxZ = std::max(mMaxZ, rhs.mMaxZ);
        }
    }IK_CONSTRAINTS;

可以看到constraints包含了yaw,pitch,roll的最大值和最小值, 作为旋转的有效范围.

在生成骨骼动画时, 提取了constraints的信息:

Vector<IK_CONSTRAINTS> constraints( boneCount );
size_t index = 0;
for(size_t i = 0; i < boneCount; ++i)
{
    scalar initPitch = 0, initYaw = 0, initRoll = 0;
    for(size_t j = 0; j < keyCountList[i]; ++j)
    {
        const BoneDQ& keyDQ = keyFrameArray[index++].getTransform();
        scalar yaw, pitch, roll;
        keyDQ.real.getYawPitchRoll(yaw, pitch, roll);
        scalar minPitch = std::min(pitch, initPitch);
        scalar maxPitch = std::max(pitch, initPitch);
        scalar minYaw = std::min(yaw, initYaw);
        scalar maxYaw = std::max(yaw, initYaw);
        scalar minRoll = std::min(roll, initRoll);
        scalar maxRoll = std::max(roll, initRoll);
        constraints[i].merge( IK_CONSTRAINTS(minPitch, maxPitch, minYaw, maxYaw, minRoll, maxRoll) );
    }
}

需要注意如果有offline工具支持skeleton文件/动画的合并操作, 那么也需要将这些约束角合并.

然后在每次CCD迭代时, 应用这些约束角. 约束角是针对关节的局部最终pose:状态量的角度, 不是单次旋转的过程量的约束.

所以在CCD中计算出的rotation, 先应用到关节点当前的pose上, 得到一个准final pose, 在用角度约束, 得到约束后的final pose:

Vector3 t = localTransformCache[i].getTranslation();
//apply rotation
Quaternion r = localTransformCache[i].getRotation() * Quaternion(axis, angle);

const IK_CONSTRAINTS& constraints = chain[i].getConstraints();
scalar yaw, pitch, roll;
//get rotated pose
r.getYawPitchRoll(yaw, pitch, roll);

//apply constraints
pitch = Math::Clamp(pitch, constraints.mMinX, constraints.mMaxX);
yaw = Math::Clamp(yaw, constraints.mMinY, constraints.mMaxY);
roll = Math::Clamp(roll, constraints.mMinZ, constraints.mMaxZ);

//final pose
localTransformCache[i].set(Quaternion(yaw, pitch, roll), t);

另外因为Blade最近把quaternion的operator* 重新定义为contatenate, 所以顺序跟以前不一样了.

奇异性问题和基于约束角的启发函数(singularity problem, and heuristic function based on constraints)

前面提到这种使用非固定转轴的方式求解CCD时会有奇异问题. 实例如下:

上图中Gizmo的位置为脚关节的目标位置, 然而对于膝关节来说, Pic和Pid已经是相同方向(都朝下), crossProduct(Pic, Pid) = vector 0 , 所以Pid是一个奇异方向(singular direction), 根据这两个向量得不到转轴. 所以腿部不会弯曲. 解决方法可以用welman的原始解法, 即使用已知固定转轴, 比如固定为x轴为转轴, 来让他弯曲(我没有尝试,原因如下:).
welman提到, 在奇异方向上, CCD的收敛速度会变慢.

而且, 即便用了welman的固定转轴的解法, 弯曲方向仍然受到约束角的限制. 比如没有约束时, 用CCD解出腿部可能能以"<"姿势达到目标, 也可能以">"姿势. 然而实际上膝盖关节不可能向外拐. 合理的解是">". 这就靠约束角来限制了, 

可是CCD在这个问题上, 在最开始迭代时, 可能就倾向于向外拐, 最后被约束, 实际上迭代中一直在尝试向外拐, 而最终结果没有转动.

比如: 把目标朝外(下图中朝右)偏移, 这个时候已经不是奇异配置了, 但是目标点偏外, 而CCD的启发方式比较激进, 是末端器离目标"越近越好", 所以迭代时, 将尝试将膝关节向外拐, 来靠近目标, 然而会被约束角限制, 其实不能转动, 只能是父节点通过类似的方式转动. 最后经过几次迭代, 又变成了奇异方向, 无法达到目标, 这种情况也要使用固定转轴:

而实际上想要的结果, 是这样 (膝关节向内)

这个时候, 如果把约束角也加入到启发过程中: 如果发现一个方向被约束, 根本行不通, 再怎么迭代下去也没意义, 那就尝试朝着不受约束角限制的方向旋转.

这个额外的启发函数, 可以避免关节朝着无意义的方向旋转. 而他恰好同时也可以将singular direction时不会旋转的问题, 变成朝着一个不会受约束角限制方向的旋转.

//fix singular direction problem. and apply heuristic direction on constraints: if impossible(angle clamped) on one direction, try the other.
if( pitch == constraints.mMinX && constraints.mMinX > -5e-2f )
    pitch = constraints.mMaxX*1e-2f;
else if( pitch == constraints.mMaxX && constraints.mMaxX < 5e-2f )
    pitch = constraints.mMinX*1e-2f;

if( yaw == constraints.mMinY && constraints.mMinY > -5e-2f )
    yaw = constraints.mMaxY*1e-2f;
else if( yaw == constraints.mMaxY && constraints.mMaxY < 5e-2f )
    yaw = constraints.mMinY*1e-2f;

if( roll == constraints.mMinZ && constraints.mMinZ > -5e-2f )
    roll = constraints.mMaxZ*1e-2f;
else if( roll == constraints.mMaxZ && constraints.mMaxZ < 5e-2f )
    roll = constraints.mMinZ*1e-2f;

这个启发值是一个非常小的偏移值, 如果没有效果, 后面迭代中会被抵消/忽略. 如果有效, 就会影响后面的迭代. 目前blade中这个启发方式比较暴力, 直接hard code, 但是原理上是这样了.

实际使用(Use IK in engine/application)

实际使用时需要设计接口, 来设置IKSover的目标. 比如通过physics引擎得到脚部的位置, 把这个位置作为腿部IK chain的目标. Blade的IK solving是在FK动画结束以后, 基于FK动画的结果来做, 所以不需要额外的blending.

关于IK chain的生成, Blade是在runtime(加载时)根据骨骼名字来建立, 而不是离线生成.

另外, Blade的IK模式也分为两种:

• Simple IK模式: 一个skeleton包含多个IK chain, 比如腿和胳膊, 4条IK chain, 但这几条IK chain是独立的, 互不影响, 也不会影响整个身体的姿势. 比如两条腿在盆骨(pelvis)处合并, 那么盆骨就作为腿部两个chain的shared base, 到这里结束. 盆骨不会参与IK计算, 所以两条腿互不影响. 这种方式十分简单, 适合用于只有手部或者脚部定位(foot placement)的需求.
  
  
• Fullbody IK模式: 整个skeleton就是一个唯一的IK chain, 它包含了多个end effector. 每个effector的定位都可能影响到身体的整个pose, 所以会影响其他effector, 比如定位左脚的时候, 盆骨也会旋转, 导致右脚受到影响. 网上很多paper里说CCD只能解决单个effector的IK chain, 对于multiple effectors, CCD不适用, 实际上我试了是可以的. 算法的整体思路是受到FABRIK(forward and backward reach inverse kinematics)的启发.(http://www.academia.edu/9165835/FABRIK_A_fast_iterative_solver_for_the_Inverse_Kinematics_problem  pdf) 当然这里实际用的是CCD解算.

Fullbody IK适合用于复杂的需求, 比如我非常喜欢的的<古墓丽影>系列, 攀爬跑跳抓, 需要用到这个方式.

  

需要注意的点: simple IK脚部定位的时候不能影响根节点, 所以如果遇到高低不平的地面, 需要寻找最低点, 把整个角色定位到最低点, 再计算IK. 而Fullbody IK理论上可以自动计算身体的位置, 但需要将跟节点设置为位移型的关节, 而不是旋转型关节, 来重新定位整个身体的位置, 不过Blade尚未支持, 后面有时间的话再加. 而且移动身体位置并不通用, 可能最为solve的参数更好, 例如: 手部reach的时候如果也移动整个身体, 可能会发生瞬移. 还有地面高度差太大, 也许要处理. 等等这些具体逻辑就不展开了.

其他遗留问题: 目前使用的是手部/脚部的定位, 而手指/脚趾的精确定位,有时候也需要, 这个以后慢慢细化.

  

误差控制 (error tolerance)

 目前Blade使用的误差是0.001, 因为IK计算是在模型空间和骨骼局部空间, 所以不需要考虑模型的缩放. 但是这个0.001是以模型/骨骼本身大小为参考, 参考值为2.0(米)高度. 对于不同大小的模型, 这个误差值会基于参考值进行缩放.

static const scalar REFERENCE_SIZESQ = (1.5f*1.5f) + (2.0f*2.0f) + (0.5f*0.5f);
static const scalar MIN_DISTANCE = 0.001f;
static const scalar MIN_DISTANCESQ = (MIN_DISTANCE*MIN_DISTANCE) / REFERENCE_SIZESQ;
const scalar tolerance = MIN_DISTANCESQ*(mIK->getSquaredSize());

实际的模型大小, 可以通过模型原始大小取得, 也可以通过骨架的binding pose的大小取得. Blade为了动画跟模型尽量解耦, 使用的是整个骨架的包围盒半径.

   

效率 (performance)

测试结果的效率还算可以, 后续可能会继续优化, 一些优化的小细节后面再记. 目前CPU i7 4770K, 单线程计算, 解算一个有效IK关节数为18, end effecotor数量为4的Full body IK, 最大迭代次数20, 最大耗时在~0.4ms左右.

而4个chain的Simple IK 解算时间大约为0.1ms.

如果配合动画LOD, 选择性的开启IK, 比如只有主角色开启IK, 或者近处角色, 可以实用.

  

记了这么多感觉有点累, 后面如果有时间再继续写备忘. IK的坑去年就打算入, 可是工作太忙. 目前mile stone 2: model算是已经结束, mile stone 3估计又到明年年底才能开始, 每年只有这段时期有点空闲. 而且工作很忙, 本身也需要休息. 后面deferred shading的计划是这样, 使用INTZ做depth pre pass, 从而充分发挥early Z的效果, 然后MRT渲染法线和颜色. dx9以后的API都可以直接读取深度缓冲, 从而不需要再单独渲染深度, 而nvidia显卡从g80以后也有INTZ来暴露新的API特性, 所以可以使用.这样既可以最大发挥earlyZ, 避免overdraw, 同时也少了G buffer的大小.

最后放一个Simple IK和Fullbody IK的对比(Fullbody时身体也会有倾斜). 另外尝试了一下用Fullbody IK定位四肢, 把stand动画改成一个攀爬动画, 也是可行的.