#include <iostream>
#include <vector>
#include <opencv2/core.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/features2d.hpp>
#include <opencv2/calib3d.hpp>
#include <opencv2/xfeatures2d.hpp>
#include <opencv2/stitching.hpp>

#define PARLIAMENT01 "x1.bmp"
#define PARLIAMENT02 "x2.bmp"

using namespace cv;
using namespace std;

int main()
{
	Mat image1 = imread(PARLIAMENT01, 0);
	Mat image2 = imread(PARLIAMENT02, 0);
	if (!image1.data || !image2.data)
		return 0;

	//namedWindow("Image 1", 1);
	//namedWindow("Image 2", 1);
	//imshow("Image 1", image1);
	//imshow("Image 2", image2);

	vector<KeyPoint> keypoints1;
	vector<KeyPoint> keypoints2;
	Mat descriptors1, descriptors2;

	//创建SIFT检测器
	Ptr<Feature2D> ptrFeature2D = xfeatures2d::SIFT::create(74);

	//检测SIFT特征并生成描述子
	ptrFeature2D->detectAndCompute(image1, noArray(), keypoints1, descriptors1);
	ptrFeature2D->detectAndCompute(image2, noArray(), keypoints2, descriptors2);

	cout << "Number of feature points (1): " << keypoints1.size() << endl;
	cout << "Number of feature points (2): " << keypoints2.size() << endl;

	//使用欧氏距离和交叉匹配策略进行图像匹配
	BFMatcher matcher(NORM_L2, true);
	vector<DMatch> matches;
	matcher.match(descriptors1, descriptors2, matches);

	Mat imageMatches;
	drawMatches(image1, keypoints1,  // 1st image and its keypoints
		image2, keypoints2,  // 2nd image and its keypoints
		matches,            // the matches
		imageMatches,       // the image produced
		Scalar(255, 255, 255),  // color of the lines
		Scalar(255, 255, 255),  // color of the keypoints
		vector<char>(),
		2);

//	namedWindow("Matches (pure rotation case)", 0);
//	imshow("Matches (pure rotation case)", imageMatches);

	//将keypoints类型转换为Point2f
	vector<Point2f> points1, points2;
	for (vector<DMatch>::const_iterator it = matches.begin();
		it != matches.end(); ++it)
	{
		float x = keypoints1[it->queryIdx].pt.x;
		float y = keypoints1[it->queryIdx].pt.y;
		points1.push_back(Point2f(x, y));

		x = keypoints2[it->trainIdx].pt.x;
		y = keypoints2[it->trainIdx].pt.y;
		points2.push_back(Point2f(x, y));
	}

	cout << "number of points: " << points1.size() << " & " << points2.size() << endl;

	//使用RANSAC算法估算单应矩阵
	vector<char> inliers;
	Mat homography = findHomography(
		points1, points2, // corresponding points
		inliers,         // outputed inliers matches
		RANSAC,      // RANSAC method
		1.);             // max distance to reprojection point

						 //画出局内匹配项
	drawMatches(image1, keypoints1,  // 1st image and its keypoints
		image2, keypoints2,  // 2nd image and its keypoints
		matches,            // the matches
		imageMatches,       // the image produced
		Scalar(255, 255, 255),  // color of the lines
		Scalar(255, 255, 255),  // color of the keypoints
		inliers,
		2);

	namedWindow("Homography inlier points", 0);
	imshow("Homography inlier points", imageMatches);

	//用单应矩阵对图像进行变换
	Mat result;
	warpPerspective(image1, // input image
		result,         // output image
		homography,     // homography
		Size(2 * image1.cols, image1.rows)); // size of output image

	cout << "homography:" << endl;
	cout << homography << endl;
											 //拼接
	Mat half(result, Rect(0, 0, image2.cols, image2.rows));
	image2.copyTo(half);
	namedWindow("Image mosaic", 0);
	imshow("Image mosaic", result);

	waitKey();
	return 0;
}

 - void mywarpPerspective(Mat src,Mat &dst,Mat T) {


   //此处注意计算模型的坐标系与Mat的不同

   //图像以左上点为（0,0），向左为x轴，向下为y轴，所以前期搜索到的特征点 存的格式是（图像x，图像y）---（rows，cols）

   //而Mat矩阵的是向下为x轴，向左为y轴，所以存的方向为（图像y，图像x）----（cols，rows）----（width，height）

   //这个是计算的时候容易弄混的


   //创建原图的四个顶点的3*4矩阵（此处我的顺序为左上，右上，左下，右下）

       Mat tmp(3, 4, CV_64FC1, 1);
       tmp.at < double >(0, 0) = 0;
       tmp.at < double >(1, 0) = 0;
       tmp.at < double >(0, 1) = src.cols;
       tmp.at < double >(1, 1) = 0;
       tmp.at < double >(0, 2) = 0;
       tmp.at < double >(1, 2) = src.rows;
       tmp.at < double >(0, 3) = src.cols;
       tmp.at < double >(1, 3) = src.rows;

   //获得原图四个顶点变换后的坐标，计算变换后的图像尺寸
       Mat corner = T * tmp;      //corner=(x,y)=(cols,rows)
       int width = 0, height = 0;
       double maxw = corner.at < double >(0, 0)/ corner.at < double >(2,0);
       double minw = corner.at < double >(0, 0)/ corner.at < double >(2,0);
       double maxh = corner.at < double >(1, 0)/ corner.at < double >(2,0);

       double minh = corner.at < double >(1, 0)/ corner.at < double >(2,0);
       for (int i = 1; i < 4; i++) {
           maxw = max(maxw, corner.at < double >(0, i)  /  corner.at < double >(2, i));
           minw  = min (minw,  corner.at < double >(0, i)  /  corner.at < double >(2, i));
           maxh  = max(maxh, corner.at < double >(1, i)  /  corner.at < double >(2, i));
           minh  = min (minh,   corner.at < double >(1, i)  /  corner.at < double >(2, i));

       }

   //创建向前映射矩阵 map_x, map_y


         //size(height,width)
       dst.create(int(maxh - minh), int(maxw - minw), src.type());
       Mat map_x(dst.size(), CV_32FC1);
       Mat map_y(dst.size(), CV_32FC1);

       Mat proj(3,1, CV_32FC1,1);

       Mat point(3,1, CV_32FC1,1);

      T.convertTo(T, CV_32FC1);  
   //本句是为了令T与point同类型（同类型才可以相乘，否则报错，也可以使用T.convertTo(T, point.type() );）
   Mat Tinv=T.inv();   

       for (int i = 0; i < dst.rows; i++) {        
           for (int j = 0; j < dst.cols; j++) {
               point.at<float>(0) = j + minw ;
               point.at<float>(1) = i + minh ;
               proj = Tinv * point;
               map_x.at<float>(i, j) = proj.at<float>(0)/ proj.at<float>(2) ;
               map_y.at<float>(i, j) = proj.at<float>(1) / proj.at<float>(2) ;
           }
       }
       remap(src,dst,map_x,map_y, CV_INTER_LINEAR);


   }