全景拼接學習-原理篇 (1) 兩張圖片之間關系計算 單應性Homograph估計

教程 https://zhuanlan.zhihu.com/p/74597564

 

 目錄
一 圖像變換與平面坐標系的關系
二 平面坐標系與齊次坐標系
三 單應性變換

 

 

一 圖像變換與平面坐標系的關系

• 旋轉：

 

 

 

 

 寫成矩陣乘法形式：

 

 

• 平移：

 

 

 

 

 但是現在遇到困難了，平移無法寫成和上面旋轉一樣的矩陣乘法形式。所以引入齊次坐標  ，再寫成矩陣形式：

 

 其中  表示單位矩陣，而  表示平移向量。

 

那么就可以把把旋轉和平移統一寫在一個矩陣乘法公式中，即剛體變換：

 

 而旋轉矩陣  是正交矩陣（  ）。

 

剛體變換：旋轉+平移（正方形-正方形）

作用：z軸距離不變 x  y 和原來相等  

仿射變換

作用：z軸距離不變 x  y 各自被比例拉伸 

 

 

 其中  可以是任意2x2矩陣（與  一定是正交矩陣不同）。

 

 仿射變換（正方形-平行四邊形）

可以看到，相比剛體變換（旋轉和平移），仿射變換除了改變目標位置，還改變目標的形狀，但是會保持物體的“平直性”。

不同  和  矩陣對應的各種基本仿射變換：

 

仿射變換（正方形-平行四邊形）

 可以看到，相比剛體變換（旋轉和平移），仿射變換除了改變目標位置，還改變目標的形狀，但是會保持物體的“平直性”。

不同  和  矩陣對應的各種基本仿射變換：

 

 

 

• 投影變換（單應性變換）

 

投影變換（單應性變換）

作用：z軸距離被拉伸   x  y 被比例拉伸  z

 

 

 

 

 投影變換（正方形-任意四邊形）

總結一下：

1. 剛體變換：平移+旋轉，只改變物體位置，不改變物體形狀  x y  z等於原來
2. 仿射變換：改變物體位置和形狀，但是保持“平直性”    z不變  x y 被比例拉伸
3. 投影變換：徹底改變物體位置和形狀                         z  x y 都被比例拉伸

 

 

 注：上圖“投影變換”應該是“任意四邊形”

我們來看看完整投影變換矩陣各個參數的物理含義：

 

 

 其中  代表仿射變換參數

  代表平移變換參數

  表示一種“變換后邊緣交點“關系，如：

 

 

至於  則是一個與  相關的縮放因子。

一般情況下都會通過歸一化使得  （原因見下文）。

二 平面坐標系與齊次坐標系

問題來了，齊次坐標到底是什么？

齊次坐標系  與常見的三維空間坐標系  不同，只有兩個自由度：

而  （其中  ）對應坐標  和  的縮放尺度。當  時：

特別的當  時，對應無窮遠：

三 單應性變換

• 單應性是什么？

此處不經證明的給出：同一個 [無鏡頭畸變] 的相機從不同位置拍攝 [同一平面物體] 的圖像之間存在單應性，可以用 [投影變換] 表示 。

注意：單應性成立是有條件的！

 

 

簡單說就是：

其中  是Left view圖片上的點，  是Right view圖片上對應的點。

• 那么這個  單應性矩陣如何求解呢？

更一般的，每一組匹配點  有

由平面坐標與齊次坐標對應關系  ，上式可以表示為：

進一步變換為：

寫成矩陣  形式：

也就是說一組匹配點  可以獲得2組方程。

 

• 單應性矩陣8自由度

注意觀察：單應性矩陣  與  其實完全一樣（其中 ），例如：

即點  無論經過  還是  映射，變化后都是  。

如果使  ，那么有:

所以單應性矩陣  雖然有9個未知數，但只有8個自由度。

在求  時一般添加約束  （也有用  約束），所以還有  共8個未知數。由於一組匹配點  對應2組方程，那么只需要  組不共線的匹配點即可求解  的唯一解。

 

 

可以看到：

1. 紅框所在平面上內容基本對齊，但受到鏡頭畸變影響無法完全對齊；
2. 平面外背景物體不符合單應性原理，偏離很大，完全無法對齊。

 

• 傳統方法估計單應性矩陣

一般傳統方法估計單應性變換矩陣，需要經過以下4個步驟：

1. 提取每張圖SIFT/SURF/FAST/ORB等特征點
2. 提取每個特征點對應的描述子
3. 通過匹配特征點描述子，找到兩張圖中匹配的特征點對（這里可能存在錯誤匹配）
4. 使用RANSAC算法剔除錯誤匹配
5. 求解方程組，計算Homograph單應性變換矩陣

樣例1  輸入兩張圖+特征檢測+自動匹配+自動求解單應矩陣+根據單應矩陣變換圖像+輸出拼接圖

　　

 opencv 349+vs2015 opencv編譯了擴展庫及其sfft庫  使用sift角點

#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;
}

　

　樣例2  變換矩陣如何實現變換  

上面是調用庫函數實現，從源碼查看自己怎么實現

dst(x,y) = src((M11x+M12y+M13)/(M31x+M32y+M33), (M21x+M22y+M23)/(M31x+M32y+M33))

 - 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);


   }