OpenCV>> camera relative pose estimation

Question

Enviroment

1. 3D Scene generated using 3DS Max
2. Camera FOV is 45 degrees
3. Two images are rendered using the same camera @ 800x600 resolution
4. Image A having camera Z rotation == 0 deg
5. Image B having camera Z rotation == 25 deg
6. 8 corresponding points ( manually drived ), no outliers

---

Task at hand

Resolve the relative camera pose between Image A & Image B ( above ) expecting to result the Induced 25 degrees rotation on the Z axis

---

Implementation

Option A:

1. Correspondence is manually generated guaranteeing no outliers ( see 'rotZ0' & 'rotZ25' at the code snippet bellow )
2. Camera Focal Length in pixel resolution is derived using Image resolution & FOV in accordance to this link
3. Camera Intrinsic matrix is composed ( according to this link ) consisting of the Image resolution and the camera FOV
4. Fundamental matrix is derived using 'cv::findFundamentalMat'
5. Essential matrix is composed ( in accordance with this link ) as a function of the Camera Intrinsic Matrix 'K' and the fundamental matrix 'F' in the following manner: 'K.t() * F * K'* where 'K.t()' is the intrinsic matrix transpose.
6. Perform SVD on the essential matrix 'matE'
7. Resolve the 4 possible solutions: [U*W*Vt], [U*W.t()*Vt], [U*W.t()*Vt.t()] & [U*W*Vt.t()]

Option B:

1. Correspondence is manually generated guaranteeing no outliers ( see 'rotZ0' & 'rotZ25' )
2. Essential Matrix is composed using 'cv::findEssentialMat'
3. Camera pose is estimated using 'cv::recoverPose'

---

Result

Neither of the above options is able to properly recover the relative camera pose (expected to be 25deg Rotation on the Z axis)

What am I doing wrong?
How can the camera relative pose be properly resolved?

Any help would be appreciated.

---

Full Code

#define RAD2DEG(rad) (((rad) * 180)/M_PI)
#define DEG2RAD(deg) (((deg) * M_PI)/180)
#define FOV2FOCAL(pixelssensorsize, fov) ((pixelssensorsize) / (2 * tan((fov) / 2)))// http://books.google.co.il/books?id=bXzAlkODwa8C&pg=PA48&lpg=PA48&dq=expressing+focal+length+in+pixels&source=bl&ots=gY4972kxAC&sig=U1BUeNHhOHmYIrDrO0YDb1DrNng&hl=en&sa=X&ei=45dLU9u9DIyv7QbN2oGIDA&ved=0CGsQ6AEwCA#v=onepage&q=expressing%20focal%20length%20in%20pixels&f=false

// http://nghiaho.com/?page_id=846
void DecomposeRotation(IN const cv::Mat& R, OUT float& fX, OUT float& fY, OUT float& fZ) {// Taken from MatLab
    fX = (float)atan2(R.at<double>(2, 1), R.at<double>(2, 2));
    fY = (float)atan2(-R.at<double>(2, 0), sqrt(R.at<double>(2, 1)*R.at<double>(2, 1) + R.at<double>(2, 2)*R.at<double>(2, 2)));
    fZ = (float)atan2(R.at<double>(1, 0), R.at<double>(0, 0));
}

int _tmain(int argc, _TCHAR* argv[])
{
    // 25 deg rotation in the Z axis (800x600)
    const cv::Point2f rotZ0[] = { { 109, 250 }, { 175, 266 }, { 204, 279 }, { 221, 253 }, { 324, 281 }, { 312, 319 }, { 328, 352 }, { 322, 365 } };
    const cv::Point2f rotZ25[] = { { 510, 234 }, { 569, 622 }, { 593, 278 }, { 616, 257 }, { 716, 303 }, { 698, 340 }, { 707, 377 }, { 697, 390 } };
    const cv::Point2f rotZminus15[] = { { 37, 260 }, { 106, 275 }, { 135, 286 }, { 152, 260 }, { 258, 284 }, { 248, 324 }, { 266, 356 }, { 260, 370 } };


    const double        dFOV = DEG2RAD(45);
    const cv::Point2d   res(800, 600);
    const cv::Point2d   pntPriciplePoint(res.x / 2, res.y / 2);
    const cv::Point2d   pntFocal(FOV2FOCAL(res.x, dFOV), FOV2FOCAL(res.y, dFOV));

    //transfer the vector of points to the appropriate opencv matrix structures
    const int                numPoints = sizeof(rotZ0) / sizeof(rotZ0[0]);
    std::vector<cv::Point2f> vecPnt1(numPoints);
    std::vector<cv::Point2f> vecPnt2(numPoints);

    for (int i = 0; i < numPoints; i++) {
        vecPnt2[i] = rotZ0[i];
        //vecPnt2[i] = rotZminus15[i];
        vecPnt1[i] = rotZ25[i];
    }

    //// Normalize points
    //for (int i = 0; i < numPoints; i++) {
    //  vecPnt1[i].x = (vecPnt1[i].x - pntPriciplePoint.x) / pntFocal.x;
    //  vecPnt1[i].y = (vecPnt1[i].y - pntPriciplePoint.y) / pntFocal.y;

    //  vecPnt2[i].x = (vecPnt2[i].x - pntPriciplePoint.x) / pntFocal.x;
    //  vecPnt2[i].y = (vecPnt2[i].y - pntPriciplePoint.y) / pntFocal.y;
    //}

    try {
        // http://docs.opencv.org/modules/calib3d/doc/camera_calibration_and_3d_reconstruction.html
        cv::Mat matK = cv::Mat::zeros(3, 3, CV_64F);
        matK.at<double>(0, 0) = pntFocal.x;
        matK.at<double>(1, 1) = pntFocal.y;
        matK.at<double>(0, 2) = pntPriciplePoint.x;
        matK.at<double>(1, 2) = pntPriciplePoint.y;
        matK.at<double>(2, 2) = 1;

        float x, y, z;
        cv::Mat R1, R2, R3, R4;
        cv::Mat t;
        cv::Mat matE;

#if 1 // Option [A]
        cv::Mat matF = cv::findFundamentalMat(vecPnt1, vecPnt2);
        matE = matK.t() * matF * matK; // http://en.wikipedia.org/wiki/Essential_matrix

        cv::Mat _tmp;
        cv::Mat U;
        cv::Mat Vt;

        cv::SVD::compute(matE, _tmp, U, Vt);

        cv::Matx33d W(0, -1, 0,
                      1,  0, 0,
                      0,  0, 1);

        R1 = U*cv::Mat(W)*Vt; // See http://stackoverflow.com/questions/14150152/extract-translation-and-rotation-from-fundamental-matrix for details
        R2 = U*cv::Mat(W)*Vt.t();
        R3 = U*cv::Mat(W).t()*Vt;
        R4 = U*cv::Mat(W).t()*Vt.t();
#else // Option [B] 
        matE = cv::findEssentialMat(vecPnt1, vecPnt2, pntFocal.x, pntPriciplePoint);// http://docs.opencv.org/trunk/modules/calib3d/doc/camera_calibration_and_3d_reconstruction.html
        cv::decomposeEssentialMat(matE, R1, R2, t);
        int iInliers = cv::recoverPose(matE, vecPnt1, vecPnt2, R4, t);// , pntFocal.x, pntPriciplePoint);
        R3 = cv::Mat::zeros(3, 3, CV_64F);
#endif

        DecomposeRotation(R1, x, y, z);
        std::cout << "Euler Angles R1 (X,Y,Z): " << RAD2DEG(x) << ", " << RAD2DEG(y) << ", " << RAD2DEG(z) << std::endl;
        DecomposeRotation(R2, x, y, z);
        std::cout << "             R2 (X,Y,Z): " << RAD2DEG(x) << ", " << RAD2DEG(y) << ", " << RAD2DEG(z) << std::endl;
        DecomposeRotation(R3, x, y, z);
        std::cout << "             R3 (X,Y,Z): " << RAD2DEG(x) << ", " << RAD2DEG(y) << ", " << RAD2DEG(z) << std::endl;
        DecomposeRotation(R4, x, y, z);
        std::cout << "             R4 (X,Y,Z): " << RAD2DEG(x) << ", " << RAD2DEG(y) << ", " << RAD2DEG(z) << std::endl;

        //cv::Mat res = matFrom.t() * matF * matTo;// Results in a null vector ( as it should ) http://en.wikipedia.org/wiki/Fundamental_matrix_(computer_vision)
        //res = matFrom.t() * matE * matTo;// Results in a null vector ( as it should )
    }
    catch (cv::Exception e) {
        _ASSERT(FALSE);
    }
    return 0;
}

---

Execution results

Option A:

Euler Angles R1 (X,Y,Z): -26.2625, 8.70029, 163.643
             R2 (X,Y,Z): 16.6929, -29.9901, -3.81642
             R3 (X,Y,Z): 5.59033, -20.841, -19.9316
             R4 (X,Y,Z): -5.76906, 7.25413, -179.086

Option B:

Euler Angles R1 (X,Y,Z): -13.8355, 3.0098, 171.451
             R2 (X,Y,Z): 2.22802, -22.3479, -11.332
             R3 (X,Y,Z): 0, -0, 0
             R4 (X,Y,Z): 2.22802, -22.3479, -11.332

Answer 1

I think that none of your results is equal to 90 deg is because of coordinate system. The coordinate system in which you rotate the camera by 90 deg about z axis is defined by yourself. You can imagine a camera pose as a vector pointing to the scene, and the R matrix decomposed from essential matrix denotes the rotation of the vector. The result is in another coordinate system. You can check my hypothesis by checking if the two camera vectors have a vector angle of 90 deg.

Answer 2

First of all, calibrate your camera instead of using predefined values. It would always make a big impact. The relative pose computed by 8-pt or 5-pt are subject to a lot of noise and by no means imply the final word. Having said that, it would be a good idea to reconstruct the points and then bundle adjust the entire scene. Refine your extrinsics and you should come up with a better pose.