圖像拼接
一、基本介紹
（image alignment）和圖像融合是圖像拼接的兩個關鍵技術。圖像拼接是計算機視覺中的重要分支，它是將兩幅以上的具有部分重疊的圖像進行無縫拼接從而得到較高分辨率或寬視角的圖像。
二、圖像拼接整體流程：
●根據給定圖像/集，實現特征匹配
●通過匹配特征計算圖像之間的變換結構利用圖像變換結構，實現圖像映射
●針對疊加后的圖像，采用APAP之類的算法，對齊特征點
● 通過圖割方法，自動選取拼接縫
●根據multi-band bleing策略實現融合
在執行RANSAC之后，我們只能在圖像中看到正確的匹配，因為RANSAC找到了一個與大多數點相關的單應矩陣，並將不正確的匹配作為異常值丟棄
（1）單應矩陣（Homography）：
有了兩組相關點，接下來就需要建立兩組點的轉換關系，也就是圖像變換關系。
單應性是兩個空間之間的映射，常用於表示同一場景的兩個圖像之間的對應關系，可以匹配大部分相關的特征點，並且能實現圖像投影，使一張圖通過投影和另一張圖實現大面積的重合。
設2個圖像的匹配點分別是$X=[x,y]^T$,$X'=[x',y']^T$，則必須滿足公式：X′=HX 且由於兩向量共線，所以X′timesHX=0
其中，$H$ 為8參數的變換矩陣，可知四點確定一個H$$begin{pmatrix}x' \y'\1 end{pmatrix} =begin{pmatrix}h{11} & h{12} & h_{13}\h{21} & h{22} & h_{23}\h{31} & h{32} & 1end{pmatrix}begin{pmatrix}x\y\1\end{pmatrix} $$
令 h=(h11:h12:h13:h21:h22:h23:h31:h32:h33)T 則有：Bh=0 N個點對給出2N個線性約束。undersethmin║Bh║2，║h║=1
（2）用RANSAC方法估算H：
1、首先檢測兩邊圖像的角點
2、在角點之間應用方差歸一化相關，收集相關性足夠高的對，形成一組候選匹配。
3、選擇四個點，計算H
4、選擇與單應性一致的配對。如果對於某些閾值:Dist(Hp、q) <ε，則點對(p, q)被認為與單應性H一致
5、重復34步，直到足夠多的點對滿足H
6、使用所有滿足條件的點對，通過公式重新計算H
7、圖像變形和融合
8、最后一步是將所有輸入圖像變形並融合到一個符合的輸出圖像中。基本上，我們可以簡單地將所有輸入的圖像變形到一個平面上，這個平面名為復合全景平面。
（3）圖像變形步驟
首先計算每個輸入圖像的變形圖像坐標范圍，得到輸出圖像大小，可以很容易地通過映射每個源圖像的四個角並且計算坐標(x,y)的最小值和最大值確定輸出圖像的大小。最后，需要計算指定參考圖像原點相對於輸出全景圖的偏移量的偏移量xoffset和偏移量yoffset。
下一步是使用上面所述的反向變形，將每個輸入圖像的像素映射到參考圖像定義的平面上，分別執行點的正向變形和反向變形。
平滑過渡（transition smoothing）圖像融合方法包括 羽化（feathering）， 金字塔（pyramid）， 梯度（gradient）
（4）圖形融合
最后一步是在重疊區域融合像素顏色，以避免接縫。最簡單的可用形式是使用羽化（feathering），它使用加權平均顏色值融合重疊的像素。我們通常使用alpha因子，通常稱為alpha通道，它在中心像素處的值為1，在與邊界像素線性遞減后變為0。
當輸出拼接圖像中至少有兩幅重疊圖像時，我們將使用如下的alpha值來計算其中一個像素處的顏色
假設兩個圖像 $I1,I2$,在輸出圖像中重疊，每個像素點$(x,y)$在圖像 I_i(x,y)=(alpha iR,alpha iG,alpha iB,alpha j,)，其中（R,G,B）是像素的顏色值，我們將在縫合后的輸出圖像中計算(x, y)的像素值：
[(α1R,α1G,α1B,α1)(α2R,α2G,α2B,α2)]/(α1α2)
三、實驗內容：
（一）圖像拼接處理：
很多錯誤匹配的特征點的，這回導致配准后的結果出現兩幅圖片對不齊的結果，所以這里要做的是刪除這些錯誤的匹配點。RANSAC(Random Sample Consensus),它是根據一組包含異常數據的樣本數據集，計算出數據的數學模型參數，得到有效樣本數據的算法。
代碼實現：
#include <iostream>
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include <opencv2/opencv.hpp>
#include<opencv2/xfeatures2d.hpp>
#include<opencv2/core/core.hpp>
using namespace cv;
using namespace std;
using namespace cv::xfeatures2d;//只有加上這句命名空間，SiftFeatureDetector and SiftFeatureExtractor才可以使用
int main()
{
//Create SIFT class pointer
Ptr<Feature2D> f2d = xfeatures2d::SIFT::create();
//SiftFeatureDetector siftDetector;
//Loading images
Mat img_1 = imread("101200.jpg");
Mat img_2 = imread("101201.jpg");
if (!img_1.data || !img_2.data)
{
cout << "Reading picture error！" << endl;
return false;
}
//Detect the keypoints
double t0 = getTickCount();//當前滴答數
vector<KeyPoint> keypoints_1, keypoints_2;
f2d->detect(img_1, keypoints_1);
f2d->detect(img_2, keypoints_2);
cout << "The keypoints number of img1 is:" << keypoints_1.size() << endl;
cout << "The keypoints number of img2 is:" << keypoints_2.size() << endl;
//Calculate descriptors (feature vectors)
Mat descriptors_1, descriptors_2;
f2d->compute(img_1, keypoints_1, descriptors_1);
f2d->compute(img_2, keypoints_2, descriptors_2);
double freq = getTickFrequency();
double tt = ((double)getTickCount() - t0) / freq;
cout << "Extract SIFT Time:" <<tt<<"ms"<< endl;
//畫關鍵點
Mat img_keypoints_1, img_keypoints_2;
drawKeypoints(img_1,keypoints_1,img_keypoints_1,Scalar::all(-1),0);
drawKeypoints(img_2, keypoints_2, img_keypoints_2, Scalar::all(-1), 0);
//imshow("img_keypoints_1",img_keypoints_1);
//imshow("img_keypoints_2",img_keypoints_2);
//Matching descriptor vector using BFMatcher
BFMatcher matcher;
vector<DMatch> matches;
matcher.match(descriptors_1, descriptors_2, matches);
cout << "The number of match:" << matches.size()<<endl;
//繪制匹配出的關鍵點
Mat img_matches;
drawMatches(img_1, keypoints_1, img_2, keypoints_2, matches, img_matches);
//imshow("Match image",img_matches);
//計算匹配結果中距離最大和距離最小值
double min_dist = matches[0].distance, max_dist = matches[0].distance;
for (int m = 0; m < matches.size(); m++)
{
if (matches[m].distance<min_dist)
{
min_dist = matches[m].distance;
}
if (matches[m].distance>max_dist)
{
max_dist = matches[m].distance;
}
}
cout << "min dist=" << min_dist << endl;
cout << "max dist=" << max_dist << endl;
//篩選出較好的匹配點
vector<DMatch> goodMatches;
for (int m = 0; m < matches.size(); m++)
{
if (matches[m].distance < 0.6*max_dist)
{
goodMatches.push_back(matches[m]);
}
}
cout << "The number of good matches:" <<goodMatches.size()<< endl;
//畫出匹配結果
Mat img_out;
//紅色連接的是匹配的特征點數，綠色連接的是未匹配的特征點數
//matchColor – Color of matches (lines and connected keypoints). If matchColor==Scalar::all(-1) , the color is generated randomly.
//singlePointColor – Color of single keypoints(circles), which means that keypoints do not have the matches.If singlePointColor == Scalar::all(-1), the color is generated randomly.
//CV_RGB(0, 255, 0)存儲順序為R-G-B,表示綠色
drawMatches(img_1, keypoints_1, img_2, keypoints_2, goodMatches, img_out, Scalar::all(-1), CV_RGB(0, 0, 255), Mat(), 2);
imshow("good Matches",img_out);
//RANSAC匹配過程
vector<DMatch> m_Matches;
m_Matches = goodMatches;
int ptCount = goodMatches.size();
if (ptCount < 100)
{
cout << "Don't find enough match points" << endl;
return 0;
}
//坐標轉換為float類型
vector <KeyPoint> RAN_KP1, RAN_KP2;
//size_t是標准C庫中定義的，應為unsigned int，在64位系統中為long unsigned int,在C++中為了適應不同的平台，增加可移植性。
for (size_t i = 0; i < m_Matches.size(); i++)
{
RAN_KP1.push_back(keypoints_1[goodMatches[i].queryIdx]);
RAN_KP2.push_back(keypoints_2[goodMatches[i].trainIdx]);
//RAN_KP1是要存儲img01中能與img02匹配的點
//goodMatches存儲了這些匹配點對的img01和img02的索引值
}
//坐標變換
vector <Point2f> p01, p02;
for (size_t i = 0; i < m_Matches.size(); i++)
{
p01.push_back(RAN_KP1[i].pt);
p02.push_back(RAN_KP2[i].pt);
}
/*vector <Point2f> img1_corners(4);
img1_corners[0] = Point(0,0);
img1_corners[1] = Point(img_1.cols,0);
img1_corners[2] = Point(img_1.cols, img_1.rows);
img1_corners[3] = Point(0, img_1.rows);
vector <Point2f> img2_corners(4);*/
////求轉換矩陣
//Mat m_homography;
//vector<uchar> m;
//m_homography = findHomography(p01, p02, RANSAC);//尋找匹配圖像
//求基礎矩陣 Fundamental,3*3的基礎矩陣
vector<uchar> RansacStatus;
Mat Fundamental = findFundamentalMat(p01, p02, RansacStatus, FM_RANSAC);
//重新定義關鍵點RR_KP和RR_matches來存儲新的關鍵點和基礎矩陣，通過RansacStatus來刪除誤匹配點
vector <KeyPoint> RR_KP1, RR_KP2;
vector <DMatch> RR_matches;
int index = 0;
for (size_t i = 0; i < m_Matches.size(); i++)
{
if (RansacStatus[i] != 0)
{
RR_KP1.push_back(RAN_KP1[i]);
RR_KP2.push_back(RAN_KP2[i]);
m_Matches[i].queryIdx = index;
m_Matches[i].trainIdx = index;
RR_matches.push_back(m_Matches[i]);
index++;
}
}
cout << "RANSAC后匹配點數" <<RR_matches.size()<< endl;
Mat img_RR_matches;
drawMatches(img_1, RR_KP1, img_2, RR_KP2, RR_matches, img_RR_matches);
imshow("After RANSAC",img_RR_matches);
//等待任意按鍵按下
waitKey(0);
}
結果顯示：
很明顯圖中的大部分錯誤點被刪除掉。接下來就可以對長焦進行矯正了，使用單映性矩陣計算，這里注意OpenCv有findFundamentalMat和findHomography兩種方法，不要搞混，為了清楚區別基本矩陣和單映性矩陣的區別，請看：
單應矩陣、基本矩陣、本質矩陣
vector<cv::Point2f> Tele_point, Wide_point;
for (int i = 0; i < InlierMatches.size(); i++)
{
Tele_point.push_back(key_points_1[InlierMatches[i].queryIdx].pt);
Wide_point.push_back(key_points_2[InlierMatches[i].trainIdx].pt);
}
cv::Mat Homography = cv::findHomography(Tele_point, Wide_point, CV_RANSAC); //計算將p2投影到p1上的單映性矩陣
cv::Mat Registration;
warpPerspective(Tele, Registration, Homography, cv::Size(Wide.cols, Wide.rows));
cv::imwrite("C:/Users/lxy/Desktop/re.jpg", Registration);
接下來將兩幅圖片融合在一起
cv::Mat Stitch(Wide.rows,Wide.cols,CV_8UC3);
Wide.copyTo(Stitch(cv::Rect(0, 0, Wide.cols, Wide.rows)));
cv::Rect Mask_center;
Mask_center.y = Wide.rows / 2 - Wide.rows / (4);//根據廣角和長焦的焦距參數決定的，lz使用的27mm和52mm的兩個鏡頭
Mask_center.x = Wide.cols / 2 - Wide.cols / (4);
Mask_center.width = Wide.cols / (2);
Mask_center.height = Wide.rows / (2);
Registration(Mask_center).copyTo(Stitch(Mask_center));
cv::imwrite("C:/Users/lxy/Desktop/Stitch.jpg", Stitch);
拼接結果：
（二）圖像拼接產生重影：
原因分析：
1.拍攝因素：在處理圖像的技術領域中，不管是對圖像進行怎樣的處理，都首先要通過拍攝設備拍取圖像，作為后續圖像處理的第一步，圖像的拍攝方式，對圖像拼接有着十分重要的影響。按照相機的拍攝情況，可以將拍攝圖像的方式分為三類。
2.圖像比例：拍攝的圖片必須要包含超過一定比例的重疊區域，當重疊區域過小時，圖像會由於匹配的特征點對不相機繞着相機的垂直軸旋轉，然后每旋轉一定的角度拍攝一張圖片。 此外拍攝的圖片必須要包含超過一定比例的重疊區域，當重疊區域過小時，圖像會由於匹配的特征點對不足，而導致拼接失敗。當重疊區域過大時，則需要多張圖像進行拼接才能獲取到更寬視野的圖像，影響拼接的實時性。一般情況下，重疊區域占整幅圖像的40%- -50%， 圖像能夠找到足夠配准的特征點對。
3.水平移動相機獲得。水平移動相機拍攝的情況是在確定要拍攝的平面后，保持相機的姿態，平行於該平面進行移動。通過這種拍攝方式拍攝出的圖片都位於同一平面上，拍攝時相機距離拍攝的目標越遠，則目標越小。由於采用水平移動相機的方式對拍攝的要求非常苛刻，因此在現實應用中很少有采用這種拍攝方進行拍攝。
小結：
理想情況下，無論何種拍攝方式，相鄰圖像的重疊區域都應該具有相同的特征，但是實際操作過程中，一方面由於拍攝圖像的存在一定的時間差，且相機的軌跡很難控制，使得圖像間光照強度可能存在差異，另一方面，拍攝圖像過程中相機存在視差，圖像重疊區域中目標並不能完全重合，此外,若拍攝場景中存在運動目標，由於拍攝時存在時間差，圖像間重疊區域中必然存在差異。在圖像拼接過程中，不管是何種原因導致了的圖像重疊區域存在不同，這必然使得最終的拼接結果存在重影現象。
（三）全景圖像拼接：
數據集：
代碼實現：
from pylab import *
from numpy import *
from PIL import Image
#If you have PCV installed, these imports should work
from PCV.geometry import homography, warp
import sift
from PCV.tools.imtools import get_imlist
"""
This is the panorama example from section 3.3.
"""
#set paths to data folder
#featname = ['C://Users//Garfield//Desktop//towelmatch//' + str(i + 1) + 'out_sift_1.txt' for i in range(5)]
#imname = ['C://Users//Garfield//Desktop//towelmatch//' + str(i + 1) + '.jpg' for i in range(5)]
imname = ['C://Users//Desktop//towelmatch//pinjie//' + str(i + 1) + '.jpg' for i in range(5)]
download_path = "C://Users//Desktop//towelmatch//pinjie//" # set this to the path where you downloaded the panoramio images
imlist = get_imlist(download_path)
l = {}
d = {}
featname = ['out_sift_1.txt' for imna in imlist]
for i, imna in enumerate(imlist):
sift.process_image(imna, featname[i])
l[i], d[i] = sift.read_features_from_file(featname[i])
#extract features and match
#l = {}
#d matches = {}
for i in range(4):
matches[i] = sift.match(d[i + 1], d[i])
#visualize the matches (Figure 3-11 in the book)
#sift匹配可視化
for i in range(4):
im1 = array(Image.open(imname[i]))
im2 = array(Image.open(imname[i + 1]))
figure()
sift.plot_matches(im2, im1, l[i + 1], l[i], matches[i], show_below=True)
#function to convert the matches to hom. points
#將匹配轉換成齊次坐標點的函數
def convert_points(j):
ndx = matches[j].nonzero()[0]
fp = homography.make_homog(l[j + 1][ndx, :2].T)
ndx2 = [int(matches[j][i]) for i in ndx]
tp = homography.make_homog(l[j][ndx2, :2].T)
# switch x and y - TODO this should move elsewhere
fp = vstack([fp[1], fp[0], fp[2]])
tp = vstack([tp[1], tp[0], tp[2]])
return fp, tp
#estimate the homographies
#估計單應性矩陣
model = homography.RansacModel()
fp, tp = convert_points(1)
H_12 = homography.H_from_ransac(fp, tp, model)[0] # im 1 to 2 # im1 到 im2 的單應性矩陣
fp, tp = convert_points(0)
H_01 = homography.H_from_ransac(fp, tp, model)[0] # = {}
tp, fp = convert_points(2) # NB: reverse order
H_32 = homography.H_from_ransac(fp, tp, model)[0] # im 3 to 2
tp, fp = convert_points(3) # NB: reverse order
H_43 = homography.H_from_ransac(fp, tp, model)[0] # im 4 to 3
#warp the images
#扭曲圖像
delta = 2000 # for padding and translation 用於填充和平移
im1 = array(Image.open(imname[1]), "uint8")
im2 = array(Image.open(imname[2]), "uint8")
im_12 = warp.panorama(H_12, im1, im2, delta, delta)
im1 = array(Image.open(imname[0]), "f")
im_02 = warp.panorama(dot(H_12, H_01), im1, im_12, delta, delta)
im1 = array(Image.open(imname[3]), "f")
im_32 = warp.panorama(H_32, im1, im_02, delta, delta)
im1 = array(Image.open(imname[4]), "f")
im_42 = warp.panorama(dot(H_32, H_43), im1, im_32, delta, 2 * delta)
imsave('C://Users//Desktop//towelmatch//result.jpg', array(im_42, "uint8"))
figure()
imshow(array(im_42, "uint8"))
axis('off')
show()
運行結果：
室內： 
室外：
全景圖像拼接選取的五張圖片由於拍攝角度的轉換存在光線問題，導致拼接不完整，不能選取比較灰暗的圖像進行。選取室外五張屋頂圖片進行運行，比最開始的五張效果好了一些，環境問題建築物少，沒有太多角特征點所以拼接效果不是很完美。
總結：
圖像拼接是把來自多個不同視角相機的圖像變換到同一視角下，無縫拼接成一張寬視野圖像（比如360度全景圖，甚至360度*180度的球面全景）。需要注意的是，由於相機各自的指向角度不一樣，因此兩圖片中來自同樣場景的部分並不能夠通過平移圖像而完全重合。事實上什么情況下不同位置、視角的相機可以變換到同一視角下拼接起來也是久為人知的，在兩種情況下圖像可拼接，一是各相機幾何中心重合；二是各相機位置任意，但場景是一個平面。一種特殊情況，場景為遠景時，可以近似的等價於一個平面場景，從而也是可拼的。在滿足上述兩種情況之一時，存在一個單應性變換，能夠將一個相機的圖像變換到另一個相機的視角下從而可以進行拼接。實際應用中為了創建出完美的全景圖，有很多的問題需要考慮。最典型的問題有兩個，一個是如何解決不同照片中曝光不一致的問題；一個是如何在拼接縫處完美平滑的融合兩張圖像的問題。第一個由曝光補償算法來解決，大體思路是估計兩張圖間的曝光差異，然后進行補償。在本次實驗過程中，不停選取圖片進行處理，對圖像特征的要求較高才能達到很好的拼接效果。
參考鏈接：https://blog.csdn.net/yangpan011/article/details/81209724