I am trying to perform space carving/voxel carving with projection matrices generated in opencv
My approach consist in:
- Camera calibration using chessboard to get the intrinsic matrix K
- 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])
- 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.
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);