This page is meant to coordinate our community efforts to implement full support for 3D Reconstruction in OpenCV.
Roadmap
Document, what's in OpenCV and what is missing
Checklist
Camera Calibration Triangulation
The simple linear method exists already in cvaux/src/cvtrifocal.cpp: icvReconstructPointsFor3View. The method can be easily adapted to two views. Someone interested in implementing the optimal algorithm 12.1 in Multiple View Geometry (2.ed) ?
Todo
Higher-level Optical Flow:
Verification via epipolar geometry reconstruction and backprojection error measurement - ransac etc.
Tri-focal tensor
Some preliminary notes
There are a lot of functions for 3D reconstuction in the cvaux folder that are not available yet but seem to work or need some small changes. For example a function for computing the projection matrices for three views and N points. From the camera matrices it is possible the define a trifocal tensor. This function uses a bundle adjustement method and functions for calculating camera matrices from point correspondences and homogeneous 4D points from camera matrices and 2D points. All this functions are not declared extern so you can not use them. To make this work you have to put the definitions in the cvaux.h. So for the mentioned function this would be:
/* Computes projection matrices for 6 points */ CVAPI(int) icvComputeProjectMatrices6Points(CvMat *points1, CvMat *points2, CvMat *points3, CvMat *projMatr1, CvMat *projMatr2, CvMat *projMatr3); /* Computes projection matrices for N points */ /* @param (input) CvMat *points1 3xN matrix for image 1 @param (input) CvMat *points2 3xN matrix for image 2 @param (input) CvMat *points3 3xN matrix for image 3 @param (output) CvMat *projMatr1 3x4 matrix for image 1 @param (output) CvMat *projMatr2 3x4 matrix for image 2 @param (output) CvMat *projMatr3 3x4 matrix for image 3 @param (input) double threshold threshold for good point @param (input) double probability probability of good result @param (output) CvMat *status 1xN matrix status of points (good or bad) @param (output) CvMat *points4D 4xN matrix of reconstructed points @result: 0 => did not find projection matrices, 1 => found projection matrices */ CVAPI(int) icvComputeProjectMatricesNPoints(CvMat *points1, CvMat *points2, CvMat *points3, CvMat *projMatr1, CvMat *projMatr2, CvMat *projMatr3, double threshold, double p, CvMat *status, CvMat *points4D);
Unfortunately there is a bug (I (StefanHafeneger) think they changed some code after they wrote this functions). To fix this you have to change three lines in cvtrifocal.cpp:
line 1571: change double result_dat[2*3]; to double result_dat[3]; line 1573: change result = cvMat(2,3,CV_64F,result_dat); to result = cvMat(1,3,CV_64F,result_dat); line 1589: change if( fabs(cvmGet(&result,1,i)) < 1e-8 ) to if( 1 == 1 )
There seem to be a lot of functions described in Multiple View Geometry in Computer Graphic by Hartley and Zisserman. So I think the work would be to debug the functions and them add the declaration in the cvaux.h.
Bugs
I (StefanHafeneger) now tried to use icvComputeProjectMatricesNPoints() for a set of tracked points. There seems to be a bug in the Bundle Adjustment. I get an assert() in cvSolveCubic(). I haven't had time yet to figure out what's going wrong there. So if someone wants to help please let me know.
Results
After applying code changes mentioned above by Stefan:
The results of icvComputeProjectMatricesNPoints() in trifocal.cpp are correct for the following synthetic test:
The test used only translation and no rotation between cameras. With this configuration the plane at infinity is determined in an easy way and hence upgrading to a metric reconstruction is possible. (Note: The first camera returned is not the canonical one. So one has to transform the cameras to apply results of the auto-calibration chapter in Multiple View Geometry.) The metric reconstruction gave little numerical error after similarity alignment to the original 3D points.
Note: There was no noise included within this test.
Further tests should include rotation between cameras and analysis of the error depending on noise.
what is in openCV
I (SebastianHoefle) start here with documentation of functions in cvcorrimages. Note, that this is primarily a documentation and I have not tested the correctness of the functions (although they look quite fine). The box at the end indicates the status of the method. So, work has to be done ...
cvaux/src/cvcorrimages.cpp
This file contains functions which track points over images and compute projection matrices for these tracked points. icvAddNewImageToPrevious adds an new image to an initial projective reconstruction and returns the projection matrix for the new image. The initial projective reconstruction can be retrieved e.g. from icvComputeProjectMatricesNPoints (in cvtrifocal.cpp) from N > 6 correspondend points in three images. (Note: In principal you could get a projective reconstruction also from two images.)
CreateFeaturePoints
This method calls cvGoodFeaturesToTrack to detect feature points.
int icvCreateFeaturePoints(IplImage *image, CvMat *points, CvMat *status)
- image
source image
- points
found feature points (type must be CV_32F or CV_64F). The desired number of features to be found is points->col.
- status
mask array for feature points (must have same number of cols than points). If the number of found features is less than points->col found, than status is zero for these points [number of foundFeatures+1, points->col]
Remarks:
should one allow to specify the parameters (e.g. minDistance, quality, blockSize ...) for cvGoodFeaturesToTrack ?
probably better to refine corner positions with cvFindCaornerSubPix. ( cvGoodFeaturesToTrack returns pixel values, not subpixels )
working
FindCorrForGivenPoints
Finds correspondend feature points with optical flow.
int icvFindCorrForGivenPoints( IplImage *image1,/* Image 1 */
IplImage *image2,/* Image 2 */
CvMat *points1,
CvMat *pntStatus1,
CvMat *points2,
CvMat *pntStatus2,
int useFilter,/*Use fundamental matrix to filter points */
double threshold);
For given points1 (with pntStatus) on image1 finds corresponding points2 on image2 and set pntStatus2 for them. The method uses optical flow ( cvCalcOpticalFlowPyrLK ) and optionally the fundamental matrix to filter the found points. The output status2 (mask mat (single channel, 8 bit)) specifies found/not found points.
Returns number of corresponding points.
working
GrowPointsAndStatus
Add to existing points and status arrays new points or just grow.
int icvGrowPointsAndStatus(CvMat **oldPoints,CvMat **oldStatus,CvMat *addPoints,CvMat *addStatus,int addCreateNum);
All points of addPoints will be added to oldPoints, in such a way that oldPoints dereferenced will point to the new allocated joined CvMat. The same for the status points. Optionally one could create more elements with argument addCreateNum (without copying, so addPoints and addStatus can be 0).
not tested
RemoveDoublePoins
Removes points which are near.
int icvRemoveDoublePoins( CvMat *oldPoints,/* Points on prev image */
CvMat *newPoints,/* New points */
CvMat *oldStatus,/* Status for old points */
CvMat *newStatus,
CvMat *origStatus,
float threshold);
Removes Points which are near ( |p - q| < threshold ) or more precisely sets origStatus to zero for the case that an old point is near to a new point. The method uses delaunay division to find nearest points between oldPoints and newPoints.
not tested
ComputeProjectMatrixStatus
Wrapper to call icvComputeProjectMatrix of cvtrifocal.cpp with the points with good status.
void icvComputeProjectMatrixStatus(CvMat *objPoints4D,CvMat *points2,CvMat *status, CvMat *projMatr)
icvComputeProjectMatrix works, and as this method copies just points with good status, this should also work.
AddNewImageToPrevious
The method adds finds new points and a projection matrix for a new image using already reconstructed 4D points.
void icvAddNewImageToPrevious____(
IplImage *newImage,//Image to add
IplImage *oldImage,//Previous image
CvMat *oldPoints,// previous 2D points on prev image (some points may be not visible)
CvMat *oldPntStatus,//Status for each point on prev image
CvMat *objPoints4D,//prev 4D points
CvMat *newPoints, //Points on new image corr for prev
CvMat *newPntStatus,// New point status for new image
CvMat *newFPoints2D1,//new feature points on prev image
CvMat *newFPoints2D2,//new feature points on new image
CvMat *newFPointsStatus,
CvMat *newProjMatr,
int useFilter,
double threshold)
This method puts things together:
first correspondend points for oldPoints are searched in newImage with icvFindCorrForGivenPoints
if more than 6 correspondences are found, the projection matrix will be computed ( calls icvComputeProjectMatrixStatus) between objPoints4D and tracked newPoints in newImage
new feature points are searched in newImage and tracked to oldImage. if the new feature points are near to the older ones found before, they are deleted with icvRemoveDoublePoins, both for the old image and the new image.
not tested
DeleteSparsInPoints
Delete Points which could not be tracked.
int icvDeleteSparsInPoints( int numImages,
CvMat **points,
CvMat **status,
CvMat *wasStatus)
This method deletes feature points, if in no image pair the correspondence could be found, i.e. the points was never found (status == 0 in all images). All points with a status of 1 in any of the images are copied to the front of the matrix and the status is set to one. The status of the remaining elements at the end is set to zero.
not tested
Outcommented methods
void icvGrowPointsArray(CvMat **points)
icvAddNewArrayPoints()
int AddImageToStruct( IplImage *newImage,//Image to add
IplImage *oldImage,//Previous image CvMat *oldPoints,// previous 2D points on prev image (some points may be not visible) CvMat *oldPntStatus,//Status for each point on prev image CvMat *objPoints4D,//prev 4D points CvMat *newPntStatus,// New point status for new image CvMat *newPoints,//New corresponding points on new image CvMat *newPoints2D1,//new points on prev image CvMat *newPoints2D2,//new points on new image CvMat *newProjMatr);//New projection matrix
These outcommented functions seem to be earlier versions of the functions above described.
New Methods
These are some new methods to support 3D reconstruction.
cv/src/cvgeometry.cpp
Compute RQ Decomposition of 3x3 matrices
/* Computes RQ decomposition for 3x3 matrices */ /* @param (input) CvMat *matrixM 3x3 matrix M @param (output) CvMat *matrixR 3x3 upper-triangular matrix R @param (output) CvMat *matrixQ 3x3 orthogonal matrix Q @param (output) CvMat *matrixQx 3x3 rotation matrix Qx (optional) @param (output) CvMat *matrixQy 3x3 rotation matrix Qy (optional) @param (output) CvMat *matrixQz 3x3 rotation matrix Qz (optional) @param (output) CvPoint3D64f *eulerAngles Euler angles of rotation (optional) */ void cvRQDecomp3x3(const CvMat *matrixM, CvMat *matrixR, CvMat *matrixQ, CvMat *matrixQx = NULL, CvMat *matrixQy = NULL, CvMat *matrixQz = NULL), CvPoint3D64f *eulerAngles CV_DEFAULT(NULL));
This function computes a RQ decomposition using Givens rotations. This function is used to decompose the left 3x3 submatrix of a projection matrix into a calibration and a rotation matrix. It optionally returns three rotation matrices, one for each axis, and the three Euler angles (for OpenGL etc.).
working
Compute projection matrix decomposition
/* Computes projection matrix decomposition */ /* @param (input) CvMat *projMatr 3x4 projection matrix @param (output) CvMat *calibMatr 3x3 internal calibration matrix @param (output) CvMat *rotMatr 3x3 external rotation matrix @param (output) CvMat *posVect 4x1 external position vector @param (output) CvMat *rotMatrX 3x3 x-axis rotation matrix (optional) @param (output) CvMat *rotMatrY 3x3 y-axis rotation matrix (optional) @param (output) CvMat *rotMatrZ 3x3 z-axis rotation matrix (optional) @param (output) CvPoint3D64f *eulerAngles Euler angles of rotation (optional) */ void cvDecomposeProjectionMatrix(const CvMat *projMatr, CvMat *calibMatr, CvMat *rotMatr, CvMat *posVect, CvMat *rotMatrX = NULL, CvMat *rotMatrY = NULL, CvMat *rotMatrZ = NULL, CvPoint3D64f *eulerAngles = NULL);
This function computes a decomposition of a projection matrix into a calibration and a rotation matrix and the position of the camera. It optionally returns three rotation matrices, one for each axis, and the three Euler angles (for OpenGL etc.).
working, implementation of infinite cameras could be added