OpenCV Stereo Camera Calibration/Image Rectification

Question

I'm trying to calibrate my two Point Grey (Blackfly) cameras for stereo vision. I'm using the tutorial stereo_calib.cpp that comes with OpenCV (code below). For some reason, I'm getting really bad results (RMS error=4.49756 and average reprojection err = 8.06533) and all my rectified images come out grey. I think my problem is that I'm not picking the right flags for the stereoCalibrate() function, but I've tried many different combinations and at best the rectified images would be warped.

Here's a link to the images I used and a sample rectified pair: https://www.dropbox.com/sh/5wp31o8xcn6vmjl/AAADAfGiaT_NyXEB3zMpcEvVa#/

Any help would be appreciated!!

static void
StereoCalib(const vector<string>& imagelist, Size boardSize, bool useCalibrated=true, bool showRectified=true)
{
    if( imagelist.size() % 2 != 0 )
    {
    cout << "Error: the image list contains odd (non-even) number of elements\n";
    return;
    }

    bool displayCorners = true;//false;//true;
    const int maxScale = 1;//2;
    const float squareSize = 1.8;
    //const float squareSize = 1.f;  // Set this to your actual square size

    // ARRAY AND VECTOR STORAGE:

    vector<vector<Point2f> > imagePoints[2];
    vector<vector<Point3f> > objectPoints;
    Size imageSize;

    //int i, j, k, nimages = (int)imagelist.size()/2;
    int i, j, k, nimages = (int)imagelist.size();

    cout << "nimages: " << nimages << "\n";

    imagePoints[0].resize(nimages);
    imagePoints[1].resize(nimages);
    vector<string> goodImageList;

    for( i = j = 0; i < nimages; i++ )
    {
    for( k = 0; k < 2; k++ )
    {
        const string& filename = imagelist[i*2+k];
        Mat img = imread(filename, 0);
        if(img.empty()) {
            break;
        }
        if( imageSize == Size() ) {
            imageSize = img.size();
        } else if( img.size() != imageSize )
        {
            cout << "The image " << filename << " has the size different from the first image size. Skipping the pair\n";
            break;
        }
        bool found = false;
        vector<Point2f>& corners = imagePoints[k][j];
        for( int scale = 1; scale <= maxScale; scale++ )
        {
            Mat timg;
            if( scale == 1 )
                timg = img;
            else
                resize(img, timg, Size(), scale, scale);
            found = findChessboardCorners(timg, boardSize, corners,
                CV_CALIB_CB_ADAPTIVE_THRESH | CV_CALIB_CB_NORMALIZE_IMAGE);
            if( found )
            {
                if( scale > 1 )
                {
                    Mat cornersMat(corners);
                    cornersMat *= 1./scale;
                }
                break;
            }
        }
        if( displayCorners )
        {
            cout << filename << endl;
            Mat cimg, cimg1;
            cvtColor(img, cimg, COLOR_GRAY2BGR);
            drawChessboardCorners(cimg, boardSize, corners, found);
            double sf = 1280./MAX(img.rows, img.cols);
            resize(cimg, cimg1, Size(), sf, sf);
            imshow("corners", cimg1);
            char c = (char)waitKey(500);
            if( c == 27 || c == 'q' || c == 'Q' ) //Allow ESC to quit
                exit(-1);
        }
        else
            putchar('.');
        if( !found ) {
        cout << "!found\n";
            break;
        }
        cornerSubPix(img, corners, Size(11,11), Size(-1,-1),
                     TermCriteria(CV_TERMCRIT_ITER+CV_TERMCRIT_EPS,
                                  30, 0.01));
    }
    if( k == 2 )
    {
        goodImageList.push_back(imagelist[i*2]);
        goodImageList.push_back(imagelist[i*2+1]);
        j++;
    }
    }
    cout << j << " pairs have been successfully detected.\n";
    nimages = j;
    if( nimages < 2 )
    {
    cout << "Error: too little pairs to run the calibration\n";
    return;
    }

    imagePoints[0].resize(nimages);
    imagePoints[1].resize(nimages);
    objectPoints.resize(nimages);

    for( i = 0; i < nimages; i++ )
    {
    for( j = 0; j < boardSize.height; j++ )
        for( k = 0; k < boardSize.width; k++ )
            objectPoints[i].push_back(Point3f(j*squareSize, k*squareSize, 0));
    }

    cout << "Running stereo calibration ...\n";

    Mat cameraMatrix[2], distCoeffs[2];
    cameraMatrix[0] = Mat::eye(3, 3, CV_64F);
    cameraMatrix[1] = Mat::eye(3, 3, CV_64F);
    Mat R, T, E, F;

    double rms = stereoCalibrate(objectPoints, imagePoints[0], imagePoints[1],
                cameraMatrix[0], distCoeffs[0],
                cameraMatrix[1], distCoeffs[1],
                imageSize, R, T, E, F,
                //TermCriteria(CV_TERMCRIT_ITER+CV_TERMCRIT_EPS, 100, 1e-5));

                TermCriteria(CV_TERMCRIT_ITER+CV_TERMCRIT_EPS, 100, 1e-5),
                CV_CALIB_FIX_ASPECT_RATIO +
                //CV_CALIB_ZERO_TANGENT_DIST +
                CV_CALIB_SAME_FOCAL_LENGTH +
                CV_CALIB_RATIONAL_MODEL +
                //CV_CALIB_FIX_K3);
                //CV_CALIB_FIX_K2);
                CV_CALIB_FIX_K3 + CV_CALIB_FIX_K4 + CV_CALIB_FIX_K5);
                //CV_CALIB_FIX_K1 + CV_CALIB_FIX_K2 + CV_CALIB_FIX_K3 + CV_CALIB_FIX_K4 + CV_CALIB_FIX_K5);
    cout << "done with RMS error=" << rms << endl;

    double err = 0;
    int npoints = 0;
    vector<Vec3f> lines[2];
    for( i = 0; i < nimages; i++ )
    {
    int npt = (int)imagePoints[0][i].size();
    Mat imgpt[2];
    for( k = 0; k < 2; k++ )
    {
        imgpt[k] = Mat(imagePoints[k][i]);
        undistortPoints(imgpt[k], imgpt[k], cameraMatrix[k], distCoeffs[k], Mat(), cameraMatrix[k]);
        computeCorrespondEpilines(imgpt[k], k+1, F, lines[k]);
    }
    for( j = 0; j < npt; j++ )
    {
        double errij = fabs(imagePoints[0][i][j].x*lines[1][j][0] +
                            imagePoints[0][i][j].y*lines[1][j][1] + lines[1][j][2]) +
                       fabs(imagePoints[1][i][j].x*lines[0][j][0] +
                            imagePoints[1][i][j].y*lines[0][j][1] + lines[0][j][2]);
        err += errij;
    }
    npoints += npt;
    }
    cout << "average reprojection err = " <<  err/npoints << endl;

    // save intrinsic parameters
    FileStorage fs("intrinsics.yml", CV_STORAGE_WRITE);
    if( fs.isOpened() )
    {
    fs << "M1" << cameraMatrix[0] << "D1" << distCoeffs[0] <<
        "M2" << cameraMatrix[1] << "D2" << distCoeffs[1];
    fs.release();
    }
    else
    cout << "Error: can not save the intrinsic parameters\n";

    Mat R1, R2, P1, P2, Q;
    Rect validRoi[2];

    stereoRectify(cameraMatrix[0], distCoeffs[0],
              cameraMatrix[1], distCoeffs[1],
              imageSize, R, T, R1, R2, P1, P2, Q,
              //CALIB_ZERO_DISPARITY, 1, imageSize, &validRoi[0], &validRoi[1]);
              CALIB_ZERO_DISPARITY, 0, imageSize, &validRoi[0], &validRoi[1]);

    fs.open("extrinsics.yml", CV_STORAGE_WRITE);
    if( fs.isOpened() )
    {
    fs << "R" << R << "T" << T << "R1" << R1 << "R2" << R2 << "P1" << P1 << "P2" << P2 << "Q" << Q;
    fs.release();
    }
    else
    cout << "Error: can not save the intrinsic parameters\n";

    // OpenCV can handle left-right
    // or up-down camera arrangements
    //bool isVerticalStereo = fabs(P2.at<double>(1, 3)) > fabs(P2.at<double>(0, 3));
    bool isVerticalStereo = false;

// COMPUTE AND DISPLAY RECTIFICATION
    if( !showRectified )
    return;

    Mat rmap[2][2];
// IF BY CALIBRATED (BOUGUET'S METHOD)
    if( useCalibrated )
    {
    // we already computed everything
    }
// OR ELSE HARTLEY'S METHOD
    else
 // use intrinsic parameters of each camera, but
 // compute the rectification transformation directly
 // from the fundamental matrix
    {
    vector<Point2f> allimgpt[2];
    for( k = 0; k < 2; k++ )
    {
        for( i = 0; i < nimages; i++ )
            std::copy(imagePoints[k][i].begin(), imagePoints[k][i].end(), back_inserter(allimgpt[k]));
    }
    F = findFundamentalMat(Mat(allimgpt[0]), Mat(allimgpt[1]), FM_8POINT, 0, 0);
    Mat H1, H2;
    stereoRectifyUncalibrated(Mat(allimgpt[0]), Mat(allimgpt[1]), F, imageSize, H1, H2, 3);

    R1 = cameraMatrix[0].inv()*H1*cameraMatrix[0];
    R2 = cameraMatrix[1].inv()*H2*cameraMatrix[1];
    P1 = cameraMatrix[0];
    P2 = cameraMatrix[1];
    }

    //Precompute maps for cv::remap()
    initUndistortRectifyMap(cameraMatrix[0], distCoeffs[0], R1, P1, imageSize, CV_16SC2, rmap[0][0], rmap[0][1]);
    initUndistortRectifyMap(cameraMatrix[1], distCoeffs[1], R2, P2, imageSize, CV_16SC2, rmap[1][0], rmap[1][1]);

    Mat canvas;
    double sf;
    int w, h;
    if( !isVerticalStereo )
    {
    sf = 600./MAX(imageSize.width, imageSize.height);
    w = cvRound(imageSize.width*sf);
    h = cvRound(imageSize.height*sf);
    canvas.create(h, w*2, CV_8UC3);
    }
    else
    {
    sf = 600./MAX(imageSize.width, imageSize.height);
    w = cvRound(imageSize.width*sf);
    h = cvRound(imageSize.height*sf);
    canvas.create(h*2, w, CV_8UC3);
    }

    for( i = 0; i < nimages; i++ )
    {
    for( k = 0; k < 2; k++ )
    {
        Mat img = imread(goodImageList[i*2+k], 0), rimg, cimg;
        remap(img, rimg, rmap[k][0], rmap[k][1], CV_INTER_LINEAR);
        cvtColor(rimg, cimg, COLOR_GRAY2BGR);
        Mat canvasPart = !isVerticalStereo ? canvas(Rect(w*k, 0, w, h)) : canvas(Rect(0, h*k, w, h));
        resize(cimg, canvasPart, canvasPart.size(), 0, 0, CV_INTER_AREA);
        if( useCalibrated )
        {
            Rect vroi(cvRound(validRoi[k].x*sf), cvRound(validRoi[k].y*sf),
                      cvRound(validRoi[k].width*sf), cvRound(validRoi[k].height*sf));
            rectangle(canvasPart, vroi, Scalar(0,0,255), 3, 8);
        }
    }

    if( !isVerticalStereo )
        for( j = 0; j < canvas.rows; j += 16 )
            line(canvas, Point(0, j), Point(canvas.cols, j), Scalar(0, 255, 0), 1, 8);
    else
        for( j = 0; j < canvas.cols; j += 16 )
            line(canvas, Point(j, 0), Point(j, canvas.rows), Scalar(0, 255, 0), 1, 8);
    imshow("rectified", canvas);
    char c = (char)waitKey();
    if( c == 27 || c == 'q' || c == 'Q' )
        break;
    }
}

Answer 1

First of all, about your calibration images. I see a few points that could lead to a better calibration :

• Use stabler images. Most of your images are blurred a bit, which results in bad accuracy in the corner detection
• Vary scale. Most images that you use present the checkerboard approx. at the same distance from the cameras.
• Be careful about your checkerboard itself. It appears to be quite badly attached to its support. If you want to achieve a good calibration, you must ensure that your checkerboard is attached tightly on a flat surface.

You have much more detailed advice about how to make a good calibration in this SO answer

Now, about the stereo calibration itself. Best way that I found to achieve a good calibration is to separately calibrate each camera intrinsics (using the calibrateCamera function) then the extrinsics (using stereoCalibrate) using the intrinsics as a guess. Have a look at the stereoCalibrate flags for how to do this.

Outside of this, your flags in the stereoCalibrate function are as such :

1. CV_CALIB_FIX_ASPECT_RATIO : you force the aspect ratio fx/fy to be fixed
2. CV_CALIB_SAME_FOCAL_LENGTH : seems OK since you have two identical cameras. You can check whether it is exact by calibrating independently each camera
3. CV_CALIB_RATIONAL_MODEL : enables K3, k4 and k5 distorsion parameters
4. CV_CALIB_FIX_K3 + CV_CALIB_FIX_K4 + CV_CALIB_FIX_K5 : fixes those 3 parameters. Since you don't use any uess, you actually put them to 0 here so the option CV_CALIB_RATIONAL_MODEL is no use in your code with those flags

Note that if you calibrate independently each camera and use the intrinsics, you have different levels of use of this data :

1. With the flag CV_CALIB_FIX_INTRINSIC, the intrinsics will be used as such and only extrinsic parameters will be optimized
2. With CV_CALIB_USE_INTRINSIC_GUESS, intrinsics will be used as guesses but optimized again
3. With a combination of CV_CALIB_FIX_PRINCIPAL_POINT, CV_CALIB_FIX_FOCAL_LENGTH and CV_CALIB_FIX_K1,...,CV_CALIB_FIX_K6 you get a little of play about which parameters are fixed and which are optimized again