C++图像拼接源程序 - - Powered by Discuz!

//根据输入图像srcImg及其角点坐标和筛选过的keypoints向量，拼接两幅图像，输出拼接后的图像及新图像的角点
class corners
{
public:
Mat left_top,left_bottom,right_top,right_bottom;
};
//图像融合总函数
void FuseImg(Mat srcImg1, corners srcImg1_in_srcImg1,
Mat srcImg2, corners srcImg2_in_srcImg2,
vector<KeyPoint> FilterdKeypoints,
vector<KeyPoint> matchedFilterdKeypoints,
Mat &dstImg, corners &dstimg_in_dstimg);
////将原图像角点变换至目标图像坐标系下，确保corners_be_transformed里的矩阵数据类型为float型
void put_value_2_transformed_corners(vector<Mat> corners_2_be_transformed,
Mat H_inverted,
vector<Mat> &corners_be_transformed);
//根据两幅图像在新图像中的位置，确定平移量
void get_transformation(vector<Mat> corner_srcImg1_in_dstImg,
vector<Mat> corner_srcImg2_in_dstImg,
Mat &trf);
//根据输入图角点在最终图中的位置向量和位移，更新输入图角点在最终图中的位置类corners，保证坐标为int型
void update_corners_in_dst(vector<Mat> corner_srcImg_in_dstImg,
Mat trf, corners &mat_src_in_dst);
//-----------------------------------------------计算dst的尺寸大小-----------------------------------------
void cal_size_dstImg(vector<Mat> corner_srcImg1_in_dstImg,
vector<Mat> corner_srcImg2_in_dstImg,
Mat trsf, int &rows_dst, int &cols_dst);
//计算像素点权值
float get_frac(int i_rows, int j_cols, corners corner_mat);
//根据输入图像角点位置更新输出图像角点
void get_union_corner(corners corner_mat_src1_in_dst, corners corner_mat_src2_in_dst, corners &corner_mat_dstImg_in_dstImg);
//将srcImg变换至dstImg上
void transf_srcImg2dst(Mat srcImg2, Mat transformation, vector<KeyPoint> FilterdKeypoints, vector<KeyPoint> matchedFilterdKeypoints, Mat &dst);
//用图像src更新图像dst
void updateDstImg(Mat srcImg,corners corner_mat_src_in_dst, corners corner_mat_dst_in_dst,Mat &dstImg);

复制代码

#include<iostream>
#include <cxcore.h>
#include <stdlib.h>
#include <time.h>
#include "opencv2/legacy/legacy.hpp"
#include <E:\xunfei\stitch_sift20200330_vc10\stitch_sift20200330_vc10\runRansac.h>
using namespace std;
using namespace cv;
/***************************** Function Prototypes ***************************/
/**
Calculates a best-fit image transform from image feature correspondences
using RANSAC.
@ keyPoints,matchedkeyPoints 特征点
@ unFilteredMatches未筛选的特征匹配
@ n 匹配数
@ m 用来计算单应性矩阵的最少特征匹配数
@param p_badxform 模型允许的最大错误概率
@param err_tol 匹配误差，通过单应性矩阵变换后，计算匹配误差，在此误差内，认为是内点
@param n_in 内点数
@filteredPoints, matchedFilteredPoints最终的内点序列
@filteredMatches，最终的特征匹配
@返回最终单应性矩阵
*/
extern void ransac_xform( vector<KeyPoint> keyPoints,
vector<KeyPoint> matchedkeyPoints,
vector<DMatch> unFilteredMatches,
int n, int m,double p_badxform,double err_tol ,
int &n_in,
vector<KeyPoint> &FilteredPoints,
vector<KeyPoint> &matchedFilteredPoints,
vector<DMatch> &FilteredMatches)
{
n=(int)unFilteredMatches.size();
Mat M;//单应性矩阵
m=N;
DMatch sample[N];//暂存被随机抽取到的匹配
DMatch* consensus;//当前一致集
DMatch* consensus_max = NULL;//当前最大一致集
double p, in_frac = RANSAC_INLIER_FRAC_EST;//p为当前计算出的模型的错误概率，当p小于p_badxform时迭代终止，in_frac为内点所占比例
//in为当前一致集元素个数，in_min为一致集中元素个数允许的最小值，保证RANSAC最终计算出的转换矩阵错误的概率小于p_badxform所需的最小内点数目
//im_max为当前最优一致集中的元素个数
//k为迭代次数
//k_max为迭代次数最大值，即初始值
int i, in, in_min, in_max = 0, k = 0, k_max;
if(n<m)
{
cout<<"实际匹配数小于计算模型需用的最小匹配数"<<endl;
}
srand( (int)time(NULL) );//初始化随机数发生器
//计算保证RANSAC最终计算出的转换矩阵错误的概率小于p_badxform所需的最小内点数目
in_min = calc_min_inliers( n, m, RANSAC_PROB_BAD_SUPP, p_badxform );
//计算k的最大值
k_max = calc_max_k(in_frac, p_badxform, m);
k=k_max;
for(int j=0;j<k;j++)
//for(int j=0;j<5;j++)
{
//pts与mpts为之后要用到的m个随机特征点对
vector<Point2f> pts,mpts;
//随机选择m(4)个匹配
draw_ransac_sample( unFilteredMatches, m, n,sample );
//提取m个匹配对应的特征点对
extract_corresp_pts( keyPoints,matchedkeyPoints,sample, m, pts, mpts );
//计算单应性矩阵M
xform_fn( pts, mpts, M );
//计算一致集和内点数
consensus = new DMatch[n];
in = find_consensus( keyPoints, matchedkeyPoints, unFilteredMatches, n, M, (double)err_tol, consensus);
//更新最优一致集的元素个数和最优一致集
if(in>in_max)
{
in_max = in;
consensus_max=consensus;
in_frac=((double)in_max)/((double)n);
}
else delete [] consensus;
//根据当前最优一致集的元素个数是否大于要求的最小内点数，判断是否继续迭代
if(in_max>in_min)
{
break;
}
//更新迭代次数k
if(calc_max_k(in_frac, p_badxform, m)>5)
{
k = calc_max_k(in_frac, p_badxform, m);
}
}
if(in_max>in_min)
{
for(int j=0; j<in_max; j++)
{
FilteredPoints.push_back(keyPoints[consensus_max[j].queryIdx]);
matchedFilteredPoints.push_back(matchedkeyPoints[consensus_max[j].trainIdx]);
FilteredMatches.push_back(consensus_max[j]);
}
n_in=in_max;
}
else
{
n_in=0;
}
}
//计算计算保证RANSAC最终计算出的转换矩阵错误的概率小于p_badxform所需的最小内点数目
static int calc_min_inliers( int n, int m, double p_badsupp, double p_badxform )
{
double pi, sum;
int i, j;
for( j = m+1; j <= n; j++ )
{
sum = 0;
for( i = j; i <= n; i++ )
{
pi = (i-m) * log( p_badsupp ) + (n-i+m) * log( 1.0 - p_badsupp ) +
log_factorial( n - m ) - log_factorial( i - m ) -
log_factorial( n - i );
/*
* Last three terms above are equivalent to log( n-m choose i-m )
*/
sum += exp( pi );
}
if( sum < p_badxform )
break;
}
return j;
}
//计算阶乘的对数
static inline double log_factorial( int n )
{
double f = 0;
int i;
for( i = 1; i <= n; i++ )
f += log((double) i );
return f;
}
//计算迭代次数k的初始值，也是最大值，程序运行后不断更新k；
int calc_max_k(double in_frac, double p_badxform, int m)
{
double p;
p=pow(in_frac, m);
p=1-p;
double s=0;
int i=0;
for(i=0;i<10000;i++)
{
s=pow(p, i);
if(s<p_badxform)
break;
}
return i;
}
//随机选择四个匹配作为计算单应性矩阵的样本值
void draw_ransac_sample( vector<DMatch> unFilteredMatches, int m, int n, DMatch* sample )
{
for(int i=0;i<m;i++)
{
int id;
for(int j=0;j<1000;j++)
{
id=rand()%n;
if(id>(((int)unFilteredMatches.size())))
break;
}
sample[i]=unFilteredMatches[id];
}
}
//提取匹配关键点KeyPoints中的点坐标输入到Point2f点对中
void extract_corresp_pts( vector<KeyPoint> kpts,
vector<KeyPoint> mkpts,
DMatch* sample, int m,
vector<Point2f> &pts,
vector<Point2f> &mpts )
{
for(int i=0;i<m;i++)
{
pts.push_back( kpts[sample[i].queryIdx].pt);
mpts.push_back( mkpts[sample[i].trainIdx].pt);
}
}
//根据4对匹配点坐标计算单应性矩阵
void xform_fn( vector<Point2f> pts, vector<Point2f> mpts, Mat &M )
{
Mat mat(8,8,CV_32FC1,Scalar::all(0));
Mat mat_invert(8,8,CV_32FC1,Scalar::all(0));
//初始化系数矩阵行0，2，4，6
for(int i_rows = 0; i_rows < 4; i_rows++)
{
mat.at<float>(2*i_rows,0) = pts[i_rows].x;//0列
mat.at<float>(2*i_rows,1) = pts[i_rows].y;//1列
mat.at<float>(2*i_rows,2) = 1.0;//2列
mat.at<float>(2*i_rows,3) = 0.0;//3列
mat.at<float>(2*i_rows,4) = 0.0;//4列
mat.at<float>(2*i_rows,5) = 0.0;//5列
mat.at<float>(2*i_rows,6) = -pts[i_rows].x*mpts[i_rows].x;//6列
mat.at<float>(2*i_rows,7) = -pts[i_rows].y*mpts[i_rows].x;//7列
}
//初始化系数矩阵行1，3，5，7
for(int i_rows = 0; i_rows<4; i_rows++)
{
mat.at<float>(1+2*i_rows,0) = 0.0;//0列
mat.at<float>(1+2*i_rows,1) = 0.0;//1列
mat.at<float>(1+2*i_rows,2) = 0.0;//2列
mat.at<float>(1+2*i_rows,3) = pts[i_rows].x;//3列
mat.at<float>(1+2*i_rows,4) = pts[i_rows].y;//4列
mat.at<float>(1+2*i_rows,5) = 1.0;//5列
mat.at<float>(1+2*i_rows,6) = -pts[i_rows].x*mpts[i_rows].y;//6列
mat.at<float>(1+2*i_rows,7) = -pts[i_rows].y*mpts[i_rows].y;//7列
}
invert(mat,mat_invert);
//计算得到向量or数组
float vct_mat_invert[8][8],vctr2[8];
for(int i_rows=0; i_rows<8; i_rows++)
{
for(int j_cols=0; j_cols<8; j_cols++)
{
vct_mat_invert[i_rows][j_cols] = mat_invert.at<float>(i_rows,j_cols);
}
}
for(int id=0; id<4;id++)
{
vctr2[id*2] = mpts[id].x;
vctr2[id*2+1] = mpts[id].y;
}
//计算homo
Mat homo(3,3,CV_32FC1,Scalar::all(0));
for(int i_rows = 0;i_rows < 3;i_rows++)
{
for(int j_cols = 0; j_cols<3 ; j_cols++)
{
if((i_rows*3+j_cols) < 8)
{
homo.at<float>(i_rows,j_cols) = calc_vct_in_prdct(vct_mat_invert[i_rows*3+j_cols], vctr2, 8);
}
else homo.at<float>(2,2)=1.0;
}
}
M=homo;
}
//计算向量（or数组）内积
float calc_vct_in_prdct(float* vctr1, float* vctr2, int dim)
{
float VectorInnerProduct=0.0;
for(int j=0;j<dim;j++)
{
VectorInnerProduct= VectorInnerProduct+vctr1[j]*vctr2[j];
}
return VectorInnerProduct;
}
//计算在单应性矩阵下，一对匹配特征点的像素距离的平方，即误差平方，返回值即为误差平方
float err_fn( KeyPoint pot, KeyPoint mpot, Mat M )
{
if((((int)M.rows)>3)||(((int)M.cols)>3))
return 1000.0;
float point[3],matched_point[3];
point[0] = pot.pt.x;
point[1] = pot.pt.y;
point[2] = 1.0;
matched_point[0] = mpot.pt.x;
matched_point[1] = mpot.pt.y;
matched_point[2] = 1.0;
float transformed_point[3] = {0.0,0.0,0.0};
for(int i=0; i<3; i++)
{
for(int j=0; j<3;j++)
{
transformed_point[i] = transformed_point[i] + M.at<float>(i,j)*point[j];
}
}
//变换后的像素坐标归一化，即将各坐标都除以Z坐标
for(int j=0; j<3; j++)
{
transformed_point[j] = transformed_point[j] / transformed_point[2];
}
float err_pixl_sqr = 0.0;
for(int i=0; i<3; i++)
{
err_pixl_sqr = err_pixl_sqr + ( matched_point[i] - transformed_point[i] ) * ( matched_point[i] - transformed_point[i] );
}
return err_pixl_sqr;
}
//寻找一致集，返回一致集元素数
int find_consensus( vector<KeyPoint> kpts, vector<KeyPoint> mkpts, vector<DMatch> uFMaches, int n, Mat M, int err_tol, DMatch* consensus)
{
int num_matches = (int)uFMaches.size();
int num_consensus=0;
double err_sqr = pow((float)err_tol, 2);
for(int i=0; i<num_matches; i++)
{
if(err_fn( kpts[uFMaches[i].queryIdx], mkpts[uFMaches[i].trainIdx], M )<err_sqr)
{
*(consensus + num_consensus)=uFMaches[i];
num_consensus++;
}
}
return num_consensus;
}

复制代码