image_rectify.cpp

基于OpenCV的计算机视觉技术实现.rar

#include "cv.h"
#include "highgui.h"

#include "cvut/cvut.h"
#include "three_dimensional_reconstruction.h"

#include <iostream>
#include <fstream>
#include <string>
#include <vector>

using namespace std;
using namespace cvut;

#pragma comment(lib,"cxcore")
#pragma comment(lib,"cv")
#pragma comment(lib,"highgui")

//#define REAL_TIME_MATCHING_SHOW

int main(int argc,char** argv)
{
	string left_image_name = "image/1_l.bmp";
	string right_image_name = "image/1_r.bmp";

	Image<uchar> left(left_image_name);
	Image<uchar> right(right_image_name);

	Matrix<double> M1(3,4,1);
	Matrix<double> M2(3,4,1);

	/************************************************************************
		获取透视投影矩阵
	*************************************************************************/
	ifstream fin("projective_matrix.txt");
	int i,j;
	for (i=0;i<3;i++) {
		for (j=0;j<4;j++) {
			fin>>M1(i,j);
		}
	}
	for (i=0;i<3;i++) {
		for (j=0;j<4;j++) {
			fin>>M2(i,j);
		}
	}

	/************************************************************************
		计算基础矩阵
	*************************************************************************/
	Matrix<double> F(calc_fundamental(M1,M2));
	
	CvMatrix3 matrix;
	for (int x=0;x<3;x++) {
		for (int y=0;y<3;y++) {
			matrix.m[x][y]=F(x,y);
		}
	}

	/************************************************************************
		图像校正
	*************************************************************************/
	cout<<"图像校正……";
	left.show("校正前 left");
	right.show("校正前 right",400,0);
	cvWaitKey();
	left.close();
	right.close();

	do_morphing(left,right,&matrix);

	left.show("校正后 left");
	right.show("校正后 right",400,0);
	cvWaitKey();
	left.close();
	right.close();
	cout<<"完成"<<endl;
	/************************************************************************
		特征提取
	*************************************************************************/
	cout<<"特征提取……";
	Image<uchar> left_gray(left.size(),8,1);
	Image<uchar> right_gray(right.size(),8,1);
	Image<uchar> left_edge(left.size(),8,1);
	Image<uchar> right_edge(right.size(),8,1);
	
	//获取灰度图
	rgb2gray(left,left_gray);
	rgb2gray(right,right_gray);

	//直方图均衡化
	hist_equalize(left_gray);
	hist_equalize(right_gray);

	//轮廓提取
 	cvCanny(left_gray.cvimage,left_edge.cvimage,10.1,10.0);
 	cvCanny(right_gray.cvimage,right_edge.cvimage,10.1,10.0);	

	left_edge.show("left");
	right_edge.show("right",400,0);
	cvWaitKey();
	left_edge.close();
	right_edge.close();
	
	//轮廓腐蚀
	cvDilate(right_edge.cvimage,right_edge.cvimage,NULL,1);
	cvDilate(left_edge.cvimage,left_edge.cvimage,NULL,1);
	left_edge.show("left");
	right_edge.show("right",400,0);
	cvWaitKey();
	left_edge.close();
	right_edge.close();

	cout<<"完成"<<endl;

	/************************************************************************
		特征匹配
	*************************************************************************/
	cout<<"特征匹配……";
	vector<vector<int> > depth;
	vector<int> depth_buf(3);
	vector<int> points_count(left.height,0);
	int win_size = 9;
	Matrix<double> left_win(1,win_size*win_size*3,1);
	Matrix<double> right_win(1,win_size*win_size*3,1);
	int depth_min = left.height;
	int depth_max = 0;
	int l,m,n;
	uchar threshold = 255;

	for (i=win_size/2;i<left.height-win_size/2;i++) {
		for (j=win_size/2;j<left.width-win_size/2;j++) {
			if (left_edge(i,j,0) >=threshold) {
				depth_buf[0] = i;
				depth_buf[1] = j;
				/* 初始最大误差 */
				double diff = 10000000.0,diff_temp;
				int index=0;
				bool is_match = false;
				for (int k=win_size/2;k<left.width-win_size/2;k++) {
					if (right_edge(i,k,0) >= threshold) {
						/* 建立原向量和目的向量 */
						int p=0,q=0;
						for (l=i-win_size/2;l<=i+win_size/2;l++) {
							for (m=j-win_size/2;m<=j+win_size/2;m++) {
								for (n=0;n<3;n++) {
									left_win(0,p++,0) = left(l,m,n);
								}
							}
							for (m=k-win_size/2;m<=k+win_size/2;m++) {
								for (n=0;n<3;n++) {
									right_win(0,q++,0) = right(l,m,n);
								}
							}
						}
						/* 计算向量误差 */
						diff_temp = cvNorm(left_win.cvmat,right_win.cvmat,CV_L2);
						/* 判断是否匹配 */
						if (diff_temp < diff){
							depth_buf[2] = k-depth_buf[1];
							diff = diff_temp;
							index = k;
							is_match = true;
						}
					}
				}
				if (is_match) {
					depth_min = abs(depth_buf[2])<depth_min?abs(depth_buf[2]):depth_min;
					depth_max = abs(depth_buf[2])>depth_max?abs(depth_buf[2]):depth_max;
					depth.push_back(depth_buf);
					is_match = false;
					points_count[i]++;
					/* 实时匹配显示 */
#ifdef REAL_TIME_MATCHING_SHOW
					cvCircle(left.cvimage,cvPoint(j,i),1,CV_RGB(255,255,255));
					cvCircle(right.cvimage,cvPoint(index,i),1,CV_RGB(255,255,255));                
					left.show("left",0,0);
					right.show("right",400,0);
					cvWaitKey();
					left.close();
					right.close();
#endif
				}
			}
		}
	}
	cout<<"完成"<<endl;
	/************************************************************************
	       显示深度图像
	*************************************************************************/
	cout<<"最小视差："<<depth_min<<endl
		<<"最大视差："<<depth_max<<endl;

	Image<uchar> depth_image(left.size(),8,1);
	for (i=0;i<depth.size();i++) {
 		depth_image(depth[i][0],depth[i][1]) = (uchar)(255*fabs((abs(depth[i][2])-depth_min)/(double)(depth_max-depth_min)));
	}

	depth_image.show("轮廓深度图");
	cvWaitKey();
	depth_image.close();

	/* 对深度图像进行线性插值 */
	int total = 0;
	for (i=0;i<points_count.size();++i) {
		for (j=0;j<points_count[i]-1;++j) {
			int step = depth_image(depth[total+j+1][0],depth[total+j+1][1])-depth_image(depth[total+j][0],depth[total+j][1]);
			step /= (depth[total+j+1][1]-depth[total+j][1]);
			for (m=depth[total+j][1]+1;m<depth[total+j+1][1];++m) {
				depth_image(depth[total+j][0],m) = depth_image(depth[total+j][0],m-1)+step;
			}
		}
		total += points_count[i];
	}

	depth_image.show("depth");
	cvWaitKey();
	depth_image.close();
//	exit(1);

	/************************************************************************
		立体显示
	*************************************************************************/
	double focal = 1000.0;
	double baseline = 2.2;
	vector<vector<double> > points3d(depth.size(),vector<double>(3));

	/* 点重建 */
	points_reconstruction(focal,baseline,depth,points3d);
	/* 三角化 */
	triangulation(points3d,depth,points_count);
	/* 三维空间显示 */
	three_dimensional_show(points3d,NULL,&argc,argv);

	return 0;
}