Kinematics describes the motion of a hierarchical skeleton structure. While forward kinematics compute world space geometric descriptions based on joint DOF values, inverse kinematics compute the vector of joint DOFs that will cause the end effector to reach some desired goal state.

In this assignment, you are required to implement two ways (Cyclic-Coordinate Descent and Jacobian Transpose) of doing inverse kinematics for an articulated character and create a demonstration that clearly shows the differences between the two techniques. We will also examine how IK works in Maya.

• Implement Cyclic-Coordinate Descent Method.
• Implement Jacobian Transpose Method.

Your handin directory is at /afs/cs.cmu.edu/academic/class/15464-s13-users/andrewid if you are in 15-464 or /afs/cs.cmu.edu/academic/class/15664-s13-users/andrewid if you are in 15-664. You should provide:

• Your code and documentation of how to run it.
• One video of your animation (two methods of IK from your code).
• A .mb file with the skeleton you load from .asf file and the iks you add to it.
• A 10-minute presentation and submit your slides or website to describe your goals, techniques you used, as well as your results.

We will use the vector: $\theta=\begin{bmatrix}\theta_1 \\ \theta_2 \\ \vdots \\ \theta_M \end{bmatrix}$ to represent the array of M joint DOF values, and the vector $e=\begin{bmatrix}e_1 \\ e_2 \\ \vdots \\ e_N\end{bmatrix}$ to represent an array of N DOFs that describes the end effector in world space. If you are using .asf skeleton file, the size of $\theta$ will depends on the joints you choose, and the size of e will be 3 since we are only concerned with the end effector position. The forward kinematics f() computes the world space end effector DOFs from the joint DOFS: $$ e = f(\theta)$$ Our goal of inverse kinematics is to compute the inverse of it: $$ \theta = f^{-1}(\theta) $$

1. Cyclic-Coordinate Descent (CCD)

CCD is a simple way to solve inverse kinematics. CCD deals with each joint individually. Although it is not as mathematically grounded as Jacobian, it's much simpler to implement. You can read these two articles "Oh My God, I Inverted Kine!" and "Making Kine More Flexible" for a thorough explanation of this method. The first article describes the general inverse kinematics concept and the second article explains the CCD method.

2. Transpose of the Jacobian

A Jacobian is a vector derivative respect to another vector. Assume $\theta$ represents the current joint DOFs and $e$ represents the current end effector DOFs, the forward kinematics function is $$e=f(\theta)$$ It maps a vector to another vector. Forward kinematics is nonlinear; it involves sin's and cos's of the input variables. Jacobian is a linear approximation of f(), it is a matrix of partial derivatives--one partial derivative for each combination of components of the vectors. We can only use the Jacobian as an approximation that is valid near the current configuration. So we must repeat the process of computing a Jacobian and taking a small step towards the goal until we get close enough. Let's use $g$ to represent the goal DOFs, here is the algorithm we can use to compute inverse kinematics:

while (e is too far from g){
    compute the Jacobian matrix $J$
    compute the transpose of the Jacobian matrix $J^T$
    compute change in joint DOFs: $\Delta \theta \approx J^T \Delta e$
    apply the change to DOFs, move a small step of $\alpha \Delta \theta$:$\theta \gets \theta + \alpha \Delta \theta$
}

Computing the Jacobian Matrix:

For any given joint pose vector $\theta$, we can explicitly compute the components of the Jacobian. $e$ is a 3D vector representing the end effector position in world space. The Jacobian is a 3xN matrix (N is the number of DOFs). Let's examine one column in the Jacobian matrix. $$ \frac{\partial e}{\partial \theta_i} = \begin{bmatrix} \frac{\partial e_x}{\partial \theta_i} & \frac{\partial e_y}{\partial \theta_i} & \frac{\partial e_z}{\partial \theta_i} \end{bmatrix}^T $$ We can fill this column numerically:
Add small $\Delta \theta_{step}$ to $\theta_i$. (Don't forget to call angles.to_pose(my_pose) after this change.) Get the change of the end effector in world space:
$\Delta e_{step} = e'-e$ where $e'$ is the new end effector position computed from angles after adding $\Delta \theta_{step}$ to $\theta_i$. (You can use the WorldBone class for this.)
So we have: $$ \frac{\partial e}{\partial \theta_i} \approx \frac{\Delta e_{step}}{\Delta \theta_{step}} = \begin{bmatrix} \frac{\Delta e_{stepx}}{\Delta \theta_{step}} & \frac{\Delta e_{stepy}}{\Delta \theta_{step}} & \frac{\Delta e_{stepz}}{\Delta \theta_{step}} \end{bmatrix}^T $$ The columns of the Jacobian matrix could also be solved analytically. See more about this in the Extra Credit section.

For Graduate Students in 15-664

For graduate students in section 15-664, "Part 1 Using the Pseudoinverse of the Jacobian Matrix" of the extra credits is mandatory.

3. IK in Maya

Set Maya back to Y-up. You can follow the instruction from this link http://mocap.cs.cmu.edu/asfamcmel/asfamcmel.php to load the skeleton from .asf file into Maya.
Click the box at the right side of "IK Handle Tool"

Change Current solver to "ikRPsolver" if it's not your default solver. (The single chain IK handle's end effector tries to reach the position and the orientation of its IK handle, whereas the rotate plane IK handle's end effector only tries to reach the position of its IK handle.)
Here's an example of how to set IK to a skeleton:
Go to side view from the viewport, select the joints of the knees, use the move tool to move them forward before we set IK. (We need to bend the knees a little bit so Maya could know which way it should bend). You can press "d" while moving, so you can move only the knee joints other than the whole hierarchy under those two joints.



You can go back to perspective view and select "IK Handle Tool", click one thigh joint, then click its corresponding ankle joint. Now you should have your IK handle set up for one leg.

Extra Credit:

1. (Mandatory for 15-664 students) Use the Pseudoinverse of the Jacobian Matrix $J^+$ instead of its transpose. Recall that the pseudoinverse of a matrix J is defined as $J^+ = J^T(JJ^T)^{-1}$. This is a closer, and more mathematically sound approximation: $$ \begin{eqnarray} \Delta e &=& J \cdot \Delta \theta \\ \Delta e &=& J I \Delta \theta \\ \Delta e &=& J (J^T) (J^T)^{-1} \Delta \theta \\ \Delta e &=& (J J^T) (J^T)^{-1} \Delta \theta \\ (J J^T)^{-1} \Delta e &=& (J^T)^{-1} \Delta \theta \\ J^T(J J^T)^{-1} \Delta e &=& \Delta \theta \\ J^+ \Delta e &=& \Delta \theta \end{eqnarray} $$ If you choose to do this method, you may want to use the safe matrix inversion code in this archive: asst2.zip
2. Interactive interface: Create a click-and-drag environment so the user can drag the target points. The character's pose updates continuously to match the targets.
   
   The best place to add interactivity is in Browser/BrowseMode.cpp. There is already an SDL event handler set up (namely "handle_event"). More thorough documentation of SDL can be found at http://www.libsdl.org/docs/html/eventfunctions.html.
3. Compute the Jacobian Matrix analytically: $$ \frac{\partial e}{\partial \theta_i} = \begin{bmatrix} \frac{\partial e}{\partial \theta_i} & \frac{\partial e}{\partial \theta_i} & \frac{\partial e}{\partial \theta_i} \end{bmatrix}^T = (a_i' \times (e-r_i')) $$ $a_i'$ is the rotation axis of this DOF in world space, $r_i'$ is the DOF's world space offset. The "axis" along with the "dof" in the asf file specify the number of joints and their axis directions. To compute the rotation axis in world space you will need this information. You need to transform them by their parent joint's world matrix.
4. Allow the user to click on any part of the character.
5. Anything else that makes posing the character more intuitive!

Tips:

Start early. This assignment is much more code and math-intensive than our 1st assignment.

Character Representations

Motion

• The heart of a Motion is the "data" field. "data" is a list of joint angles for each frame of an animation.
• The first N values are the joint angles for frame 0, the next N values are frame 1, etc. (Here N is the number of Degrees of Freedom in the skeleton)

The base code mainly treats motions as read-only, but if you plan to edit motions a lot, you may find this useful:

void Motion::set_angles(unsigned int frame, Character::Angles angles)
{
    assert(frame < frames());
    assert(skeleton == angles.skeleton);
    for(unsigned int i=0; i<data.size(), i++)
    {
        data[frame * skeleton->frame_size + i] = angles.angles[i];
    }
}
    

Angles

Angles is a rather tricky representation. It is effectively a list of joint angles, with one caveat: the first 3 values are actually the root x y and z position. After that, the next 3 values are the XYZ root rotation, followed by the XYZ rotation of bone 1, etc.

root_x = angles[0]; //position
root_y = angles[1]; //position
root_z = angles[2]; //position
root_X = angles[3]; //rotation
root_Y = angles[4]; //rotation
root_Z = angles[5]; //rotation
left_hip_X = angles[6]; //local rotation
left_hip_Y = angles[7]; //local rotation
left_hip_Z = angles[8]; //local rotation
...

This is the same representation Motion uses for each frame.

Pose

Pose is a more intuitive and higher level representation, consisting of: a root position, a root orientation, and a list of orientations (quaternions) for the other bones. To find the orientation of a particular bone (other than root joint), you'll need to know its index is in the list. (This goes for WorldBones as well). The list of bones and their corresponding id's are as follows:

lhipjoint (id=0)
lfemur (id=1)
ltibia (id=2)
lfoot (id=3)
ltoes (id=4)
lowerback (id=5)
upperback (id=6)
thorax (id=7)
lclavicle (id=8)
lhumerus (id=9)
lradius (id=10)
lwrist (id=11)
lhand (id=12)
lfingers (id=13)
lthumb (id=14)
lowerneck (id=15)
upperneck (id=16)
head (id=17)
rclavicle (id=18)
rhumerus (id=19)
rradius (id=20)
rwrist (id=21)
rhand (id=22)
rfingers (id=23)
rthumb (id=24)
rhipjoint (id=25)
rfemur (id=26)
rtibia (id=27)
rfoot (id=28)
rtoes (id=29)

For the names of the joints, see the following image.

Quatf quat;
quat = my_pose.bone_orientations[16]; //Get orientation of left shoulder
//Compute new rotation.
...
my_pose.bone_orientations[16] = quat; //Set new rotation for left shoulder

World Bones

• World bones gives you a list of joint positions at the base and the tip, and a list of orientations (quaternion).
• Build a WorldBones from a Pose
• There is no function to change WorldBones to Poses or Angles, so treat this representation as "read only" Hint: getting the tip position from WorldBones is very useful for IKs!

Character::WorldBones world_bones;
get_world_bones(my_pose, world_bones);
Vector3f tip = world_bones.tips[22]; //World position of left finger tip

Math

• 3d rotations represented by quaternions
• Make a quaternion by calling "rotation"
      -between 2 vectors
      -a rotation axis and an angle (in radians)
• Apply a rotation by calling "rotate"
• Useful properties to know:
      -multiplying 2 quaternions composes the rotations one after the other
      -conjugation changes the direction of rotation (so -q works as expected)
      -quaternion multiplication is NOT commutative

Quick reference

The base code has a lot of files to search through. Here is a list of the files and functions you're likely to find useful.

Library/Library.hpp
void Motion::get_angles(unsigned int frame, Character::Angles &into) const;
void Motion::get_pose(unsigned int frame, Character::Pose &into) const;

Character/Character.hpp
void Pose::to_angles(Angles &into) const;
void Angles::to_pose(Pose &into) const;

Character/pose_utils.hpp
void get_world_bones(Pose const &pose, WorldBones &out);