2020년 11월 17일 화요일

Inverse Kinematics(IK) Iterate Method 정리

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.97.816&rep=rep1&type=pdf 

https://medium.com/unity3danimation/overview-of-jacobian-ik-a33939639ab2

http://graphics.cs.cmu.edu/nsp/course/15-464/Spring11/handouts/iksurvey.pdf

참고할자료


1. 배경


손의 위치는 팔꿈치, 어깨 등의 관절의 각도와 각 뼈의 길이에 영향을 받는데, 원하는 위치에 손을 갖다대려면 각 관절이 어떻게 움직여야 하는가?

IK 는 위와 같은 문제를 푸는데 사용된다.

풀이 방법은 크게 두개로 수치해석으로 푸는 방법과 컴퓨터를 이용한 반복을 통해서 푸는 방법이 있다. 

전자는 수식 몇개로 끝나서 빠른데 Bone 이 많으면 수식을 확장하기가 힘들고,

후자는 확장성은 좋은데 반복 알고리즘 특으로 정확도와 속도랑 반비례한다.

후자는 CCD(Cyclic Coordinate Descent) 과 Jacobian 을 많이 쓰는데 위 논문은 Jacobian 쓰는 논문이다.


CCD 랑 Jacobian 의 큰 차이점이라면 CCD 는 각 IKChain(각 관절들) 이 한번에 하나씩만 변환되고, Jacobian 은 모든 IKChain 이 한꺼번에 변환된다는 것을 꼽을 수 있을 것이다.

한번에 하나씩면 변환되므로 CCD 는 안정성이 떨어지고 Jacobian은 안정적이다.

하지만 속도를 생각하면...



2. CCD


1. 설명


MMD(MikuMiku Dance) 에서 사용하고 언리얼, 유니티에서도 사용하는 방법이다. 

휴리스틱이라 정말 빠른데 결과물이 종종 자연스럽지 못하다.

관절에 각도 제한을 두어서 이를 해결하는데, 그래도 불완전한 경우가 꽤 많다.



아이디어는 위 그림처럼 간단하다.

ikNode 와 knee, foot 과 knee 사이의 각도를 구해서 knee 를 회전시킨다.

똑같이 hip joint 도 회전시키고, 그리고 이걸 계속 반복하면 foot 은 ikNode 에 있게 된다.


2. 구현

코드를 보면서 자세히 살펴보자.

크게 3단계로 나눴다.


1
2
3
4
5
6
7
8
for (auto& chain : ikChains)
{
    chain.prevAngle = 0;
    chain.rot = { 0,0,0,1 };
    auto ik_node = chain.node;
    ik_node->SetIKRotation(chain.rot);
    ik_node->Update();
}
cs

ikChains 를 초기화하는 단계이다. 

매번 자세를 바꿀때 마다 시행한다. 

IK 문제가 답이 여러개인 경우가 많아서, 초기 자세가 똑같아야지 같은 target 위치에 같은 결과물이 나오기 때문이다.


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
for (int i = 0; i < ikProblem.iterateCount; i++)
{
    for (auto& chain : ikProblem.ikChains)
    {
        Matrix invChain = chain.node->GetWorldMatrix().Inverse_SRT();
 
        // chainNode 의 위치가 원점인 bone space 로 ikNode(원하는 위치), targetNode(현재 위치)를 변환함
        auto ikpos = Matrix::MatrixToPosition(ikNode->GetWorldMatrix());
        auto targetpos = Matrix::MatrixToPosition(targetNode->GetWorldMatrix());
 
        // chain node 가 원점이므로, 원점과 ikNode, 원점과 TargetNode 의 벡터를 구함 
        Vector3 ikVec = (invChain * ikpos).Normalize();
        Vector3 targetVec = (invChain * targetpos).Normalize();
 
        // 위에서 구한 두 벡터같의 각도를 구함
        float theta = std::acos(Math::Clamp(Vector3::Dot(ikVec, targetVec), -1.f, 1.f));
 
        // theta 가 작으면 할일이 없음
        if (theta < 1.0e-5f)
            continue;
        
        RotateChain(theta, ikVec, targetVec);
    }
}
cs

ikNode(원하는 위치) 와 IKChain 의 노드,  TargetNode(현재 위치) 와 IKChain 의 노드 사이의 벡터를 구한다.

그리고 ikChain 을 localSpace 상에서 회전시킬 것이기 때문에 편의상 ikChain 의 BoneSpace 로 변환한다.

그렇게해서 구한 ikVec, targetVec 을 unit vector 로 만들어서 내적후 arccos 에 집어넣으면 우리가 찾고 있는 각도가 나온다. 

tan2 를 써서 하는 방법 등 구하는 법은 여러가지니 참고.

이 각도를 이용해 회전하고 이를 각 ikChain 마다 순차적으로 시행하고 이걸 또 반복한다.


이때 ikNode 와 TargetNode 의 길이가 줄어드는지 확인해서 안줄어들면 break 하는 등 무의미한 반복을 회피해서 속도를 올릴 수 있다.


1
2
3
4
5
6
7
8
9
10
void RotateChain(float theta, Vector3 ikVec, Vector3 targetVec)
{
    auto rotateAxis = Vector3::Cross(targetVec, ikVec).Normalize();
    // cross 는 순서가 있어서 순서 바꾸면 안댐
 
    auto rot = chain.rot * Quaternion::QuaternionFromAngleAxis(theta, rotateAxis);
 
    chain.node->SetIKRotation(rot);
    chain.node->Update();
}
cs

그리고 마지막 회전부분이다.

회전 축은 위에서 구한 ikVec 과 targetVec 의 cross 를 통해서 구할 수 있다.

이 축을 이용해 위에서 구한 theta 랑 같이 quaternion 을 만들어 회전시키면 된다.

 

만약 회전축이 고정되어 있다면 상황에 따라 다르지만 X, Y, Z 축에 대한 고정이면 간단하다.

회전방향을 알 수 없기 때문에 +theta, -theta 를 모두 구해서 원하는 결과물에 가까운걸 택하면 된다. 

 


3. Jacobian(Half)


1. Jacobian 이 어떻게 사용되는가


Jacobian Matrix 는 비선형 변환을 극소구간에는 선형변환으로도 가능하다는 아이디어이다.

행렬의 미분과도 비슷한데, 이걸 이용해서 gradient descent 처럼 반복해 우리가 원하는 각도를 구하는 방법이 이 방법이다.

자세히 살펴보자.


현재 IK 노드들의 각도 상태를 s 라고 하고 (ex 30, 60, 90) s 의 원하는 결과값을 r 라고 하면

r = s + t

라고 표현할 수 있다. 

t 는 위 식 그대로 s 가 r 이 되기위해서 필요한 변화량이다. 


t 를 구하는 방법은 현재 End Node 의 위치인 E와 원하는 위치 T 를 통해서 알 수 있다.

(T-E) = M * t

각도가 t 만큼 변환하면 End Node 는 (T-E) 만큼 이동할 것이고 즉 원하는 위치에 도달할 것이라는 식이다.

이때 각도가 거리로 변환되는 것은 M 이 비선형변환이기 때문에 가능한 일이다.

비선형변환이라는게 쉽게 말하면 M 에는 t 를 인자로 받아드리는 삼각함수 등이 들어있을테니 우리가 아는 행렬의 규칙들을 쓸 수 없다는 말이다.

참고로 동차행렬도 비선형변환이다.

이때 M 이 정확하게 무엇이냐하고 하면 그건 잘 모르겠고, 여기선 있다고 하고 넘어가자.



왼쪽은 각 노드의 각인 θ에 대한 공간, 오른쪽은 Global Space 에 대한 공간으로, 간단하게 2차원으로 표현했지만 실제론 다차원일 것이다.

왼쪽 공간의 격자가 비선형변환인 M 혹은 Jacobian 변환행렬들의 모음을 통해서 오른쪽 공간의 구불구불한 이상한 격자로 변환된걸 보자.

그림처럼 비선형변환은 복잡하다.


하지만 Jacobian 행렬을 사용하면 비선형대신 선형변환을 사용할 수 있다.

위에서 빨간색 선의 마디는 dθ(왼쪽), 그리고 그것에 Jacobian 을 사용한 선형변환의 결과물(오른쪽)이다.

이때 시작하는 구간은 원점이고 계속 결과물(빨간 마디)을 더해서 원하는 목적지에 다가가는 것을 확인할 수 있다.


덧붙여 선형변환이므로 빨간 선의 조각들은 직선이다.

오른쪽 공간에서 빨간 선을 구불하게 그렸지만, 선의 사이가 매우 작다면 (T-E) 를 향한 직선에 가까워질 것이다.



위 그림을 식으로 표현하면

$\left(T-E\right)=\int _{\theta =0}^{\theta }J\cdot d\theta $$\left(T-E\right)=\int _{\theta =0}^{\theta }\combi{J}_θ\cdot d\theta $(TE)=θθ=0Jθ·dθ

이라고 대략 나타낼 수 있다.


이를 간략히

dp = J * dθ 

라고 표현하면 중간에 노드가 A, B, C 가 있다고 할 때

라고 표현할 수 있다.

내부의 미분관련은 chain rule 을 이용한 Jacobian 의 식이 응용된 것이다. 

원래 공간의 델타를 변환할 공간의 델타로 바꾸고 있음을 알 수 있다.


그런데 어차피 Global Space 에선 point 는 ( T - E ) 벡터로 정해져있다.

E 는 매번 반복할 때 마다 업데이트 되어야하겠지만 말이다. 

그러므로 p_x, p_y, p_z 는 ( T - E ) 로 바꿀 수 있다.


dp_x / dθ_A 에 대해선, 정확한 증명은 공부를 안해서 모르겠으나

각 Node 에서의 rotate axis 를 이용해서 구할 수 있다.

Rotate Axis 와 ( E - Pos_ChainNode ) 의 cross 를 구하면 각 공간 축인 x, y, z 에 대한 변화율이 나온다.

머리속으로 그려보면 대충 말이 되는거 같은데...


그럼 이제 dθ 는 쉽게 구할 수 있다.

J 는 선형변환이니 역행렬을 구해서 양쪽에 곱하면 되기 때문이다.

이때 정확도를 조금 포기해서 Transpose 를 사용하면 매우 빠르다.

Psuedo Inverse 나 DLS 방법 등 여러방법이 있는데 복잡해서 여기선 생략한다.

$\begin{bmatrix}T-E\end{bmatrix}\ =\ \begin{bmatrix}R_A\ \times \left(T-E\right)&R_B\ \times \left(T-E\right)&R_C\ \times \left(T-E\right)\end{bmatrix}\ \begin{bmatrix}\combi{\ d\theta }_A\\\combi{\ d\theta }_B\\\combi{\ d\theta }_C\end{bmatrix}$[
TE
]
 = [
RA ×(TA)RB ×(TB)RC ×(TC)
]
 [
 dθA
 dθB
 dθC
]
$\ \begin{bmatrix}R_A\ \times \left(T-E\right)&R_B\ \times \left(T-E\right)&R_C\ \times \left(T-E\right)\end{bmatrix}^T\cdot \begin{bmatrix}T-E\end{bmatrix}\ =\ \combi{\ d\theta }$ [
RA ×(TA)RB ×(TB)RC ×(T− C)
]
T·[
TE
]
 = 
 dθ

Transpose(J) * ( T - E ) 를 통해 dθ 구해 ikchains 를 회전시키는 것을 반복한다.

그럼 언젠가 T - E 가 epsilon 보다 작아질 것이다.


2. 구현


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
for (int i = 0; i < ikProblem.iterateCount ; i++)
{
    for (auto& chain : ikProblem.ikChains)
    {        
        auto chainNode = _transform_array[chain.index];
        Vector3 chain_pos = chainNode->GetPosition();    
        
        for (int i = 0; i < 3; i++)
        {
            Vector3 axis = GetRotateAxis(i);
            Vector3 local_axis = Matrix::QuaternionToMatrix(chainNode->GetRotation()).Transpose() * axis;
 
            auto cross = Vector3::Cross(axis, (E - chain_pos));
            float theta = Vector3::Dot(cross, (T - E));
            chain.rot = chain.rot * Quaternion::QuaternionFromAngleAxis(theta * _delta, local_axis);
        }    
    }
 
    for (auto& chain : ikProblem.ikChains)
    {
        chain.node->SetIKRotation(chain.rot);
        chain.node->Update();
    }
}
cs

위는 jacobian 의 역행렬을 transpose 로 대체한 것이다.

theta 를 구할 때 Dot 을 한 것은 Jacobian 을 transpose 하면 cross 가 행이 되기 때문이다.

그렇게 구한 결과값을 axis 를 chain node 의 local 공간으로 변환해서 적용한다.

이때 delta 는 상황에 따라 적절한걸 택해야한다.


이때 DOF 가 1이 아닌 chain node 인 경우 축을 여러개로 분할해서 사용하면 된다.

예를들어 DOF 가 3 이면 X, Y, Z 축 에 대해서 하나의 Node 에 모두 적용하면 된다.

이는 위에서 GetRotateAxis() 로 구현했는데, i 가 3인 이유는 X, Y, Z 가 3개이기 때문이다.

List