1

I'm trying to apply IK to my human skeleton's arm, where the shoulder position is fixed, and only the arm is moving.

The FK of the arm looks like this

(x,y,z) = Root_T * Shoulder_T * Shoulder_Rx(q_0) * Shoulder_Ry(q_1) * 
Shoulder_Rz(q_2) * Elbow_T * Elbow_Rx(q_3) * Hand_T * Hand_Rx(q_4) * Position

This is my IK solver function

//Q_desired is a 5*1 Vector whose entries are initialized as 0.0
//TargetP is where I want the end effector to be
void BoneRig::solveIK(glm::vec3& targetP, Eigen::VectorXf& Q_desired) {
    //Basic setups
    float error = 10.0f;
    glm::vec3 errorP;
    Eigen::VectorXf deltaQ(5);
    glm::vec3 endEffector = glm::vec3(0.0f,0.0f,0.0f);
    Eigen::MatrixXf J(3,5);
    Eigen::MatrixXf J_plus(5,3);
    float deltaTime = 0.1;

    int limit = 1000;
    int limitCounter = 0;

    while(error > 0.01){
        //gets the current endEffector location and stores it in endEffector
        getEndEffectorP(endEffector);
        errorP = targetP - endEffector;

        Eigen::VectorXf deltaX(3);
        deltaX << (errorP * deltaTime).x, (errorP * deltaTime).y, (errorP * deltaTime).z;

        getJacobian(J, targetP);
        J_plus = J.transpose();
        deltaQ = J_plus * deltaX;

        Q_desired += deltaQ;
        setOrientation(Q_desired);

        error = glm::length(errorP);

        if(limitCounter>= limit){
            std::cout << "You could not solve this"<<std::endl;
            break;
        }
    }
}

I think the problem is in setOrientation(Q_desired); function. Q_desired = Q_desired + deltaQ updates the DOF angles in every loop, so I update the rotation matrices of each joints accordingly. Here is the setOrientation(Q_desired) function.

//Joints[15] is the Shoulder joint, Joints[16] is the Elbow joint.
//Joints[17] is the Hand joint
void BoneRig::setOrientation(const Eigen::VectorXf& ceta_d){
    glm::mat4 temp = glm::mat4(1.0f);
    Joints[15].setR(glm::rotate(temp,ceta_d(0),glm::vec3(1.0f,0.0f,0.0f)));
    Joints[15].setR(glm::rotate(Joints[15].returnR(),ceta_d(1),glm::vec3(0.0f,1.0f,0.0f)));
    Joints[15].setR(glm::rotate(Joints[15].returnR(),ceta_d(2),glm::vec3(0.0f,0.0f,1.0f)));
    temp = glm::mat4(1.0f);
    Joints[16].setR(glm::rotate(temp,ceta_d(3),glm::vec3(1.0f, 0.0f, 0.0f)));
    temp = glm::mat4(1.0f);
    Joints[17].setR(glm::rotate(temp,ceta_d(4),glm::vec3(1.0f, 0.0f, 0.0f)));
    return;
}

I was curious if applying the angle rotations should be done in those rotation axis (Should it be the global axis, and not the local axis?) I didn't post the getJacobian() method because I think it's right, but I will if there's no flaw in this code. Any help would be appreciated. Thanks in advance!


EDIT This is how I get the Jacobian matrix for the arm. Note that only the shoulder is the ball and socket joint and the rest are just revolute joints based on the local x-axis.

void BoneRig::getJacobian(Eigen::MatrixXf& J, glm::vec3& targetWorld){
    glm::vec4 xAxis = glm::vec4(1.0f,0.0f,0.0f,0.0f);
    glm::vec4 yAxis = glm::vec4(0.0f,1.0f,0.0f,0.0f);
    glm::vec4 zAxis = glm::vec4(0.0f,0.0f,1.0f,0.0f);
    glm::mat4 temp = Joints[0].returnT() * Joints[0].returnR() * Joints[11].returnT() * Joints[11].returnR() * Joints[14].returnT() * Joints[14].returnR() * Joints[15].returnT() * Joints[15].returnR();
    glm::vec3 worldXAxis = glm::vec3(temp * xAxis);
    glm::vec3 worldYAxis = glm::vec3(temp * yAxis);
    glm::vec3 worldZAxis = glm::vec3(temp * zAxis);
    glm::vec3 worldPos = glm::vec3(temp * glm::vec4(0.0f,0.0f,0.0f, 1.0f));
    glm::vec3 p = targetWorld - worldPos;

    glm::vec3 upperPart = glm::cross(worldXAxis,p);
    J(0,0) = upperPart.x;
    J(1,0) = upperPart.y;
    J(2,0) = upperPart.z;

    upperPart = glm::cross(worldYAxis,p);
    J(0,1) = upperPart.x;
    J(1,1) = upperPart.y;
    J(2,1) = upperPart.z;

    upperPart = glm::cross(worldZAxis,p);
    J(0,2) = upperPart.x;
    J(1,2) = upperPart.y;
    J(2,2) = upperPart.z;

    temp = temp * Joints[16].returnT() * Joints[16].returnR();
    worldXAxis = glm::vec3(temp * xAxis);
    worldPos = glm::vec3(temp * glm::vec4(0.0f,0.0f,0.0f, 1.0f));
    p = targetWorld - worldPos;
    upperPart = glm::cross(worldXAxis,p);
    J(0,3) = upperPart.x;
    J(1,3) = upperPart.y;
    J(2,3) = upperPart.z;
    temp = temp * Joints[17].returnT() * Joints[17].returnR();

    worldXAxis = glm::vec3(temp * xAxis);

    worldPos = glm::vec3(temp * glm::vec4(0.0f,0.0f,0.0f, 1.0f));

    p = targetWorld - worldPos;
    upperPart = glm::cross(worldXAxis,p);

    J(0,4) = upperPart.x;
    J(1,4) = upperPart.y;
    J(2,4) = upperPart.z;
}

I used this equation for each column of the Jacobian Matrix. The axis should be calculated in global coordinates, so that's what I did in the code.

enter image description here

genpfault
  • 51,148
  • 11
  • 85
  • 139
Peter
  • 460
  • 6
  • 23
  • 2
    The rotation matrices should be as defined in your forward model. Hence, this looks correct. Why do you use the transpose as the inverse of the Jacobian? This does not look correct. – Nico Schertler Oct 27 '18 at 12:31
  • @Nico Schertler But I heard that Jacobian transpose can be used as a pseudo inverse of a Jacobian, and it will still give you the answer... Is it wrong? – Peter Oct 27 '18 at 13:52
  • Only in very special cases. This does not seem to be one of them. – Nico Schertler Oct 27 '18 at 15:21

1 Answers1

0

Mr. Schertler was right. Everything in the code was correct, except that I should not have used J_plus = J.transpose. Instead, I should've used the Moore-Penrose Pseudoinverse which is

J_plus = (J.transpose() * J).inverse() * J.transpose() // This is for Overdetermined systems
J_plus = J.transpose * (J*J.transpose()).inverse() // This is for Underdetermined systems

In this case, since it's a 5DOF system with only 3 rows of the Jacobian matrix, I should've used the first version.

Thanks again for @Nico Schertler!

Peter
  • 460
  • 6
  • 23