Projection matrix for voxel carving

edit

I am trying to perform space carving/voxel carving with projection matrices generated in opencv

My approach consist in:

1. Camera calibration using chessboard to get the intrinsic matrix K
2. I use solvePNP to obtain rotation and translation vectors, rodrigues to get the rotation matrix. I asssamble R and T into one matrix which I multiply by intrinsic matrix (K*[R|T])
3. In opencv I draw a cube on the chessboard using the drawContours function and I use the projectPoints function so the cube rotates and translates together with the chessboard This seems to work- see the attached image 1.

4. Than I perform voxel carving:

   cube of 200x200x200 voxels of size 4 is generated

   foreach voxel I calculate its image reprojection [u,v,1]T= K[R|T][voxel.x,voxel.y,voxel.z,1]T

   I mask all the voxels that are ~ white (background)

After one image processing (one reprojection) I see a piramide - see figure 2.

I would expect a cube. I assumed that since projectPoints works the projection matrix is OK but maybe I am wrong or I should modify it some how before I spent days trying to figure it out with no success so if somebody could help me or point a direction that would be great

My intrinsic matrix : [5.4133176190360371e+002, 0., 3.1423029943981362e+002; 0., 5.4314440666700341e+002, 2.4109085795979527e+002; 0., 0., 1.]

First projection matrix: [-130.3838689234658, 543.2310276120221, -282.2779745141828, 4963.465064332162; 378.6239093571447, -7.628680906981948, -457.9480175480329, 4526.52325286986; -0.4457741775970285, 0.006323269925829211, -0.8951231193780763, 22.58023681026416]

generation of projection matrix

            int board_w = 9;
    int board_h = 6;
    int board_n = board_w*board_h;
    vector<Point2f> corners1;
    vector<Point3f> obj;
    Mat rvec, tvec;
    for (int j = 0; j < board_n; j++)
    {
        obj.push_back(Point3f(j / board_w, (j%board_w), 0.0));         
    }
    solvePnP(obj, corners1, intrinsics, distortion, rvec, tvec, false);
    Rodrigues(rvec, dst);
            hconcat(dst, tvec, proj);
    P = intrinsics*proj;

reprojection of the voxel cube (the code is from http://nghiaho.com/?p=2124)

        X2.at<float>(0) = vox.x;
        X2.at<float>(1) = vox.y;
        X2.at<float>(2) = vox.z;
        X2.at<float>(3) = 1.0;

        X = P * X2;

        int x = X.at<float>(0) / X.at<float>(2);
        int y = X.at<float>(1) / X.at<float>(2);