承接前文骨骼動(dòng)畫(huà)理論及程序部分實(shí)現(xiàn)（一）前向動(dòng)力學(xué)。
簡(jiǎn)介
??如果說(shuō)前向動(dòng)力學(xué)是找到父骨骼的空間，來(lái)變換得到子骨骼空間，那么通過(guò)子骨骼得到父骨骼的位置，就稱(chēng)之為反向動(dòng)力學(xué)，也稱(chēng)為逆向動(dòng)力學(xué)。
??常見(jiàn)應(yīng)用于機(jī)器人，或者即時(shí)演算游戲中，人物腿部在凹凸不平地面中的表現(xiàn)。諸如刺客信條等擁有豐富攀爬動(dòng)畫(huà)系統(tǒng)的游戲?qū)K會(huì)有更多的應(yīng)用。
IK解算
基本原理
??對(duì)于只有兩個(gè)骨骼關(guān)節(jié)點(diǎn)，其中一點(diǎn)固定，那么剩下一點(diǎn)的位置也隨之確定：
??對(duì)于有三個(gè)關(guān)節(jié)點(diǎn)的問(wèn)題，可以用余弦定理解決
??常見(jiàn)的IK骨，例如腿，用這個(gè)就能解決，不過(guò)有時(shí)候給頭發(fā)或者蜈蚣、觸手、機(jī)械臂等物體加IK時(shí)，IK鏈會(huì)很長(zhǎng)，余弦定理就無(wú)法很好的解決了。
??對(duì)于更通用的IK解算方法，常用的有循環(huán)坐標(biāo)下降法（Cyclic Coordinate Descent, 簡(jiǎn)稱(chēng)CCD）和雅克比矩陣法(Jacobian Matrix)。
??針對(duì)于MikuMikudance的模型文件存儲(chǔ)參數(shù)，很容易發(fā)現(xiàn)它是用CCD來(lái)解決IK解算問(wèn)題的，因此我們也用CCD來(lái)進(jìn)行IK解算。
循環(huán)坐標(biāo)下降法
名詞
??首先我們確定一些名詞。
- IK鏈（Link\Chain）：IK解算的單位，受IK效應(yīng)器影響的一系列父子骨骼。
- 開(kāi)始節(jié)點(diǎn)（Start Joint）：IK鏈固定的那一端節(jié)點(diǎn)。
- 末端效應(yīng)器（End Effector）：IK鏈的末端節(jié)點(diǎn)，也是盡量要到達(dá)目標(biāo)的節(jié)點(diǎn)。
- IK目標(biāo)（Target）：希望末端效應(yīng)器到達(dá)的目標(biāo)點(diǎn)。
??注意：這里我們約定，IK鏈條從接近末端開(kāi)始計(jì)數(shù)，如下圖。
??現(xiàn)在比照MMD模型對(duì)應(yīng)上面的名詞：
??其中IK鏈包含左膝（左ひざ）、開(kāi)始節(jié)點(diǎn)左足（左腿跟），以及末端效應(yīng)器。
??注意末端效應(yīng)器在MMD中其實(shí)是上面標(biāo)著Target的左足首，而真正的IK目標(biāo)其實(shí)是左足IK本身。
??如果看骨骼親子繼承關(guān)系，左足>左ひざ>左足首才是一條鏈的，而真·IK目標(biāo)的繼承關(guān)系是全ての親>左足IK親>左足IK，這意為著左足IK直屬于全親骨，只隨著全親骨的變化而變化，其位置不會(huì)輕易改變。
算法
??從IK鏈靠近末端效應(yīng)器的第一個(gè)節(jié)點(diǎn)開(kāi)始，算當(dāng)前節(jié)點(diǎn)到末端效應(yīng)器的向量與當(dāng)前節(jié)點(diǎn)到IK目標(biāo)向量的夾角，然后旋轉(zhuǎn)之；一直算到開(kāi)始節(jié)點(diǎn)，然后再?gòu)牡谝粋€(gè)節(jié)點(diǎn)開(kāi)始，不斷循環(huán)。
??如上圖，先像紅色線條那樣算角度，然后讓末端效應(yīng)器落在當(dāng)前節(jié)點(diǎn)和IK目標(biāo)的向量（下一張圖的黑色線）之上。（由于我作圖時(shí)是目測(cè)旋轉(zhuǎn)，后兩張圖沒(méi)有準(zhǔn)確落到黑色線條張，這個(gè)別在意）
??然后沒(méi)有然后了，原理就這么簡(jiǎn)單，下面開(kāi)始放部分程序?qū)崿F(xiàn)。
程序?qū)崿F(xiàn)
??首先，因?yàn)椴粩嘁鹿趋赖男D(zhuǎn)，之前前向動(dòng)力學(xué)的BNode不夠用，我們此時(shí)不光是要記錄生成骨骼空間變換本身，同時(shí)要計(jì)算它自身的旋轉(zhuǎn)：
struct BNode {
BNode(PMX::Bone& _bone) : bone(_bone) {};
int32_t index;
PMX::Bone& bone;
BNode* parent = nullptr;
std::vector<BNode*> childs;
bool haveAppendParent = false;
BNode* appendParent = nullptr;
float appendWeight;
std::vector<BNode*> appendChilds;
//棄用Mconv矩陣，而是每次用getLocalMat生成
glm::vec3 position;
glm::quat rotate;
//記錄骨骼本身經(jīng)歷的移動(dòng)和旋轉(zhuǎn)
glm::vec3 translate;
glm::quat quaternion;
inline glm::mat4 getLocalMat() {//conv before local to after world 變換初始狀態(tài)骨骼空間坐標(biāo)到完成位置的世界坐標(biāo)
return glm::translate(glm::mat4(1), position) * glm::mat4_cast(rotate);
}
inline glm::mat4 getGlobalMat() {//conv before world to after world 變換初始狀態(tài)世界坐標(biāo)到完成位置的世界坐標(biāo)
return glm::translate(glm::mat4(1), position) * glm::mat4_cast(rotate) * glm::translate(glm::mat4(1), -bone.position);
}
//更新rotate 和 position ，使用前確保所有父骨骼都是最新更新的
//只管自己，不管子骨和賦予親
bool updateSelfData() {
glm::mat4 parentLocalMat(1);
glm::vec3 parentPosition(0);
glm::quat parentRotate(1, 0, 0, 0);
glm::vec3 appendParentTranslate(0);
glm::quat appendParentQuaternion(1, 0, 0, 0);
if (parent != nullptr) {
parentPosition = parent->bone.position;
parentRotate = parent->rotate;
parentLocalMat = parent->getLocalMat();
}
if (appendParent != nullptr) {
appendParentTranslate = appendParent->translate;
appendParentQuaternion = appendParent->quaternion;
}
position = parentLocalMat * glm::vec4(bone.position - parentPosition + translate + appendParentTranslate * appendWeight, 1);
rotate = parentRotate * quaternion * glm::quat(glm::eulerAngles(appendParentQuaternion) * appendWeight);
return true;
}
//更新自身和所有子骨和賦予親，使用前確保父骨都是更新過(guò)的
//更新規(guī)則：先更新自身，然后更新自身的賦予親本身和賦予親的直屬子骨，再更新自身子骨
bool updateSelfAndChildData() {
std::stack<Animation::BNode*> nodes;
nodes.push(this);
while (!nodes.empty()) {
Animation::BNode* curr = nodes.top();
nodes.pop();
curr->updateSelfData();//更新自身
for (Animation::BNode* appendChild : curr->appendChilds) {//更新賦予親
std::stack<Animation::BNode*> appendChildNodes;
appendChildNodes.push(appendChild);
while (!appendChildNodes.empty()) {
Animation::BNode* appendChildCurr = appendChildNodes.top();
appendChildNodes.pop();
appendChildCurr->updateSelfData();
for (auto child = appendChildCurr->childs.rbegin(); child != appendChildCurr->childs.rend(); ++child) {
appendChildNodes.push(*child);
}
}
}
for (auto child = curr->childs.rbegin(); child != curr->childs.rend(); ++child) {//所有孩子壓棧
nodes.push(*child);
}
}
return true;
}
}
??吐個(gè)槽，之前沒(méi)有存translate和quaternion，每次都用親骨反向計(jì)算，后來(lái)還有顧及賦予親的影響。更新也是將當(dāng)前骨骼到所有子骨、賦予子骨都更新一遍，為了保證更新的子骨用的是親骨更新前的變換矩陣，還搞了個(gè)雙棧，分別前序和后續(xù)遍歷，結(jié)果輸出程序結(jié)果不正確，我也要被繁雜的邏輯搞炸了。完全搞不清到底是IK解算的問(wèn)題，還是更新骨骼的問(wèn)題。
??然后多存了那兩個(gè)變量，更新骨骼簡(jiǎn)化為只更新自身，程序就變得特別清爽和可控；說(shuō)多了都是淚。
??關(guān)于updateSelfAndChildData的更新規(guī)則也是根據(jù)情況來(lái)的。我發(fā)現(xiàn)賦予親和賦予子骨的父骨一般是一樣的，大致情況如下：
??結(jié)果個(gè)別骨骼就特殊了，它只繼承被賦予骨，不繼承普通骨骼，這樣的更新規(guī)則更新不到它，于是我只能設(shè)定更新賦予骨的同時(shí)，向下更新一級(jí)直接子骨。
??雖說(shuō)這樣很依賴(lài)模型本身，但也是沒(méi)辦法的事情，假如兩個(gè)骨骼互相認(rèn)爹（反向?qū)嬍谊P(guān)系），那一開(kāi)始個(gè)骨骼樹(shù)解析也是沒(méi)完沒(méi)了的。建模不規(guī)范，程序兩行淚啊。
??前向動(dòng)力學(xué)部分遍歷骨骼樹(shù)，要記錄一開(kāi)始的translate和quaternion這個(gè)沒(méi)有技術(shù)含量，代碼我就不放了（別忘了賦予親就行），接下來(lái)是重點(diǎn)的IK解算部分。
??我將IK解算器封裝為一個(gè)類(lèi)，在此之前先定義一個(gè)內(nèi)部結(jié)構(gòu)IKChain(這個(gè)命名可能有問(wèn)題，應(yīng)該叫IKJoint更準(zhǔn)確)
struct IKChain {
Animation::BNode* node;
bool enableAxisLimit;
glm::vec3 min;
glm::vec3 max;
IKChain* pre = nullptr;//前一條鏈節(jié)
Animation::BNode* endPre = nullptr;//末端效應(yīng)器沒(méi)有鏈節(jié)的屬性，額外存儲(chǔ)
bool updateChain() {//更新IK鏈
IKChain* curr = this;
do {
curr->node->updateSelfData();//更新自身
if (curr->pre == nullptr) {//如果是第一個(gè)鏈節(jié)，就更新末端效應(yīng)器
curr->endPre->updateSelfData();
}
curr = curr->pre;
} while (curr != nullptr);
return true;
}
};
??注意：這是更新順序，不是解算順序，可以理解為：第一次解算更新前一個(gè)關(guān)節(jié)就行，第二次要更新前兩個(gè)……
class IKSolve
{
public:
IKSolve(BNode* _ikNode, BoneManager& boneManager);//初始化
~IKSolve();
void Solve();//解算
private:
Animation::BNode* ikNode;
std::vector<IKChain> ikChains;
Animation::BNode* targetNode;
float limitAngle;
};
IKSolve::IKSolve(Animation::BNode* ikNode, Animation::BoneManager& boneManager)
{
assert(ikNode->bone.isIK());//確保是IK骨
targetNode = boneManager.linearList[ikNode->bone.ikTargetBoneIndex];//按照MMD命名為T(mén)arget，其實(shí)是末端效應(yīng)器
limitAngle = ikNode->bone.ikLimit;//單位角，不讓一次旋轉(zhuǎn)角度過(guò)大的限制
this->ikNode = ikNode;//IK目標(biāo)節(jié)點(diǎn)
ikChains.resize(ikNode->bone.ikLinks.size());
for (int i = 0; i < ikChains.size(); ++i) {
ikChains[i].node = boneManager.linearList[ikNode->bone.ikLinks[i].boneIndex];
ikChains[i].enableAxisLimit = ikNode->bone.ikLinks[i].enableLimit();
ikChains[i].min = ikNode->bone.ikLinks[i].min;
ikChains[i].max = ikNode->bone.ikLinks[i].max;
if (i != 0) {
ikChains[i].pre = &ikChains[i - 1];
}
}
ikChains[0].endPre = targetNode;
}
??接下來(lái)是重頭戲：
void IKSolve::Solve() {
for (int ite = 0; ite < ikNode->bone.ikIterationCount; ++ite) {
for (size_t chainIdx = 0; chainIdx < ikChains.size(); ++chainIdx) {
IKChain& chain = ikChains[chainIdx];
glm::mat4 worldToLocalMat = glm::inverse(chain.node->getLocalMat());
glm::vec3 jointLocalIK = worldToLocalMat * glm::vec4(ikNode->position, 1); //關(guān)節(jié)本地坐標(biāo)系下IK目標(biāo)位置
glm::vec3 jointLocalTarget = worldToLocalMat * glm::vec4(targetNode->position, 1); //關(guān)節(jié)本地坐標(biāo)系下末端效應(yīng)器位置
if (glm::distance(jointLocalIK, jointLocalTarget) < 1e-3) {//兩個(gè)向量差距太小，目的達(dá)到，直接退出解算
ite = ikNode->bone.ikIterationCount;
break;
}
float cos_deltaAngle = glm::dot(glm::normalize(jointLocalIK), glm::normalize(jointLocalTarget));
if (cos_deltaAngle > 1 - 1e-6) {//夾角太小pass
continue;
}
float deltaAngle = glm::acos(cos_deltaAngle);
deltaAngle = glm::clamp(deltaAngle, -limitAngle, limitAngle);//一次旋轉(zhuǎn)的度數(shù)不得超過(guò)限制角
glm::vec3 axis = glm::normalize(glm::cross(glm::normalize(jointLocalTarget), glm::normalize(jointLocalIK)));//旋轉(zhuǎn)軸
chain.node->quaternion *= glm::quat_cast(glm::rotate(glm::mat4(1), deltaAngle, axis));
chain.updateChain();
}
}
ikChains.end()[-1].node->updateSelfAndChildData();//從IK鏈末端（開(kāi)始節(jié)點(diǎn)）更新所有子骨
}
??第一重循環(huán)的循環(huán)次數(shù)模型自己會(huì)給出，第二重循環(huán)是對(duì)每個(gè)關(guān)節(jié)鏈進(jìn)行對(duì)比、旋轉(zhuǎn)。
??關(guān)于那個(gè)逆矩陣是這樣的：我們上次推導(dǎo)了矩陣變換，其中
是點(diǎn)所處骨骼空間的坐標(biāo)，
是點(diǎn)所處世界空間的坐標(biāo)，左右兩邊的左側(cè)都乘上
的逆矩陣得到新的等式：
，也就是可以輸入世界坐標(biāo)，得到骨骼關(guān)節(jié)空間的坐標(biāo)，這樣旋轉(zhuǎn)的角度能保證是基于骨骼關(guān)節(jié)坐標(biāo)系的。
??兩個(gè)向量、角度檢測(cè)不得不做，當(dāng)夾角太小時(shí)，cos值接近1，求arccos可能會(huì)出現(xiàn)nan。
??是不是很簡(jiǎn)單？（才怪，這程序?qū)懙念^都禿了）把IK解算放到FK之后（上一章POSE的構(gòu)造函數(shù)最后）：
for (Animation::BNode* node : boneManager.linearList) {
if (node->bone.isIK()) {
Animation::IKSolve solve(node, boneManager);
solve.Solve();
}
}
??左側(cè)由剛剛理論寫(xiě)成的程序生成，末端效應(yīng)器看來(lái)是落在目標(biāo)點(diǎn)上了，不過(guò)小姐姐你的左腿是不是有點(diǎn)不對(duì)勁啊，趕緊去砍醫(yī)生吧。
??才怪了，單單是我們上面那個(gè)三節(jié)點(diǎn)模型，用余弦定理，在三維空間中也能解除無(wú)數(shù)個(gè)解（第二個(gè)關(guān)節(jié)點(diǎn)在空間中一個(gè)圓邊上都可以），類(lèi)似腿這樣的關(guān)節(jié)點(diǎn)都是有限制的，膝蓋關(guān)節(jié)只能在自己關(guān)節(jié)點(diǎn)空間的X軸旋轉(zhuǎn)，并且只能向內(nèi)彎折。
??看其他人的程序，以前的PMD模型文件可能給骨骼一個(gè)標(biāo)記，告訴程序：這個(gè)骨骼是膝蓋（isLeg），也有部分程序直接把左右膝蓋的shift-jis編碼值寫(xiě)死在程序里，特殊處理這些關(guān)節(jié)。
??不過(guò)PMX文件的IK鏈中，給了每個(gè)鏈節(jié)的限制范圍（在最上面的圖能看到，注意下上面記錄的都是左手坐標(biāo)系的范圍），由此，我們對(duì)這樣的關(guān)節(jié)特殊處理：
if (chain.enableAxisLimit) {//如果這個(gè)關(guān)節(jié)點(diǎn)有軸限制
glm::mat4 inv = glm::inverse(chain.node->parent->getLocalMat());
glm::vec3 selfRotate = glm::eulerAngles(chain.node->quaternion);//本身基于父空間的旋轉(zhuǎn)的歐拉角表示
if (ite == 0) {//第一次迭代，直接到旋轉(zhuǎn)到可容忍區(qū)間的一半
glm::vec3 targetAngles = (chain.min + chain.max) / 2.0f;//目標(biāo)旋轉(zhuǎn)
chain.node->quaternion = glm::quat(targetAngles);//令關(guān)節(jié)和全部子骨旋轉(zhuǎn)
chain.updateChain();
}
else {//不是第一次迭代，也要保證旋轉(zhuǎn)在區(qū)間內(nèi)
glm::vec3 axis = glm::normalize(glm::cross(glm::normalize(jointLocalTarget), glm::normalize(jointLocalIK)));
glm::vec3 deltaRotate = glm::eulerAngles(glm::quat_cast(glm::rotate(glm::mat4(1), deltaAngle, axis)));
deltaRotate = glm::clamp(deltaRotate, chain.min - selfRotate, chain.max - selfRotate);
chain.node->quaternion *= glm::quat(deltaRotate);
chain.updateChain();
}
continue;
}
else {//沒(méi)有軸限制，程序和上邊一樣
glm::vec3 axis = glm::normalize(glm::cross(glm::normalize(jointLocalTarget), glm::normalize(jointLocalIK)));
chain.node->quaternion *= glm::quat_cast(glm::rotate(glm::mat4(1), deltaAngle, axis));
chain.updateChain();
}
??第一次迭代中，直接將關(guān)節(jié)旋轉(zhuǎn)到允許空間的一半，如膝蓋的限制是[(0~180), (0, 0), (0,0)]，那么就直接旋轉(zhuǎn)到[90, 0, 0](右手坐標(biāo)系的)。隨后其他迭代中，保證限制關(guān)節(jié)不出允許區(qū)間，例如本身是[30, 0, 0]，那么它只能旋轉(zhuǎn)范圍為：[(-30, 150), (0, 0), (0, 0)]。
其他
??IK解算參考了很多github上現(xiàn)有的代碼，但還是寫(xiě)的頭快禿了，賦予親中間還過(guò)來(lái)?yè)v亂，簡(jiǎn)直……給你們看看我中間調(diào)試時(shí)出現(xiàn)得到N種奇葩狀態(tài)：
??模型來(lái)自女仆麗塔-洛絲薇瑟2.0-神帝宇，希望不要打我（滑稽。
??如果還是有些看不懂，可以看看以下的網(wǎng)站：
引用
循環(huán)坐標(biāo)下降（CCD）算法中對(duì)骨骼動(dòng)畫(huà)中膝蓋等關(guān)節(jié)的特殊處理
挺久前一個(gè)先輩的文章，也是搞MMD的
saba-OpenGL Viewer (OBJ PMD PMX)
從我未接觸圖形學(xué)時(shí)就看上了這個(gè)github項(xiàng)目，現(xiàn)在越寫(xiě)越心驚，骨骼動(dòng)畫(huà)、表情動(dòng)畫(huà)、物理都有，以后也要繼續(xù)借鑒這里的代碼
MikuMikuDance PMX/VMD Viewer for Windows, OSX, and Linux
上面的saba項(xiàng)目太龐大，中間封裝了很多層，疑惑的時(shí)候看看輕量些的代碼也是個(gè)不錯(cuò)的選擇
MMDAgent
也是個(gè)不錯(cuò)的借鑒
