圖片重疊區(qū)域識別有很多應(yīng)用場景，例如全景照片的合成等。這篇文章介紹一種圖片重疊區(qū)域識別的方法。
1. 整體流程
識別重疊區(qū)的步驟包含：
- 查找兩張圖片的特征點；
- 匹配特征點；
- 計算單應(yīng)矩陣；
- 基于單應(yīng)矩陣，計算第二張圖片的四個頂點在圖片1中的位置；
- 計算 (4) 中得到的位置與圖片邊緣的交點，得到重疊區(qū)域的頂點；
- 同理得到圖片1的頂點在圖片2中覆蓋區(qū)域的頂點。
來看代碼和效果：
2. 查找和匹配特征點
找特征點的方式有很多，這里用的是 SIFT，也可以用 SURF、ORB 都行。目的是為了計算出單應(yīng)矩陣。
Matcher 選的是 FlannBasedMatcher， 選擇 BruteForceMatcher 也是可以的。
int main() {
// 讀取圖片
auto photo1 = cv::imread("../IMG_1841.jpeg");
auto photo2 = cv::imread("../IMG_1846.jpeg");
// 分別檢測兩張圖片的特征點和描述子
cv::Mat descriptor1, descriptor2;
std::vector<cv::KeyPoint> keyPoint1, keyPoint2;
cv::Ptr<cv::SiftFeatureDetector> detector = cv::SiftFeatureDetector::create(300);
detector->detectAndCompute(photo1, cv::Mat(), keyPoint1, descriptor1);
detector->detectAndCompute(photo2, cv::Mat(), keyPoint2, descriptor2);
// 特征點匹配
cv::FlannBasedMatcher matcher;
std::vector<cv::DMatch> matches;
matcher.match(descriptor1, descriptor2, matches);
// ① ↓
經(jīng)過這兩步我們就能得到匹配的特征點了，但這些匹配對有很多不正確的：
匹配后的特征點
通過漢明距離過濾掉部分匹配對：
// ① ↓
// 獲取匹配對漢明距離的最小值和最大值
double minDistance = 1000, maxDistance = 0;
for (auto& match : matches) {
double distance = match.distance;
if (distance < minDistance) { minDistance = distance; }
if (distance > maxDistance) { maxDistance = distance; }
}
// 對匹配對按照漢明距離篩選
std::vector<cv::DMatch> goodMatches;
for (auto& match : matches) {
if (match.distance <= 5 * minDistance) {
goodMatches.push_back(match);
}
}
// ② ↓
過濾后的匹配對：
挑選后的匹配對
不同的照片這一步得到的匹配對數(shù)量不一樣，通常變化不大的場景留下的匹配對會更多。
3. 單應(yīng)矩陣和重疊區(qū)域的計算
單應(yīng)矩陣約束了同一3D空間點在兩個像素平面的2D齊次坐標。有了單應(yīng)矩陣我們就能計算圖片2的頂點在圖片1中應(yīng)該在什么位置。
// ② ↓
// 找出關(guān)鍵點對應(yīng)的像素坐標
std::vector<cv::Point2f> points1, points2;
for (auto &match : goodMatches) {
points1.push_back(keyPoint2[match.trainIdx].pt);
points2.push_back(keyPoint1[match.queryIdx].pt);
}
// 計算單應(yīng)矩陣
cv::Mat H = findHomography(points1, points2, cv::RANSAC);
// 拿到圖片的四個頂點
std::vector<cv::Point> vertex = {
{0, 0},
{0, photo2.rows},
{photo2.cols, photo2.rows},
{photo2.cols, 0}
};
// 用于保存 照片2 的頂點 在 照片1 中的位置
std::vector<cv::Point> vertices2AtPhoto1;
// 遍歷 圖片2 的四個頂點，根據(jù)單應(yīng)矩陣，計算出這四個頂點在 圖片1 中的位置
for (auto& point : pointsB) {
cv::Point temp;
// 這個方法在下面??
transformPoint(H, point, temp);
vertices2AtPhoto1.push_back(temp);
}
// ④ ↓
其中 transformPoint 為：
/**
* 計算坐標 in 經(jīng)過 單應(yīng)矩陣 homography 變換回去的坐標，并保存在 out 中。
*/
void transformPoint(const cv::Mat& homography, const cv::Point& in, cv::Point& out) {
// 頂點的齊次坐標
cv::Vec3d homogeneous = {(double)in.x, (double)in.y, 1.0};
// 計算這個點在圖1上的位置
auto temp = (cv::Mat)(homography * homogeneous);
cv::Vec3d point = temp.col(0);
auto x = point[0];
auto y = point[1];
auto z = point[2];
out.x = cvRound(x / z);
out.y = cvRound(y / z);
}
這樣我們就得到了圖2的頂點在圖1中的位置，但是這個位置可能已經(jīng)超出圖片1的視野外了，所以想要計算重疊區(qū)域，還需要計算圖2變換后的頂點與圖一的交點。交點圍成的多邊形就是重疊區(qū)域。這是個純數(shù)學問題：
// ④ ↓
// 計算 圖像B 的頂點在圖像A中 和 圖像A邊框 的交點
std::vector<cv::Point> outIntersection1;
getPolyIntersections(pointsB, vertices2AtPhoto1, outIntersection1);
// ⑤ ↓
其中 getPolyIntersections 及相關(guān) callee 方法我放在了文章最后。
同樣，我們計算出圖1在圖2中的多邊形頂點：
// ⑤ ↓
// 用同樣的方式，計算出圖片1的四個頂點，在圖片2中與圖片2的交點。
auto hInverses = H.inv();
std::vector<cv::Point> outIntersection2;
// 遍歷 圖片2 的四個頂點，根據(jù)單應(yīng)矩陣，計算出這四個頂點在 圖片1 中的位置
for (auto& point : vertex) {
cv::Point temp;
transformPoint(hInverses, point, temp);
outIntersection2.push_back(temp);
}
// ⑥ ↓
為了更好的顯示和后續(xù)覆蓋度的計算等，我們可以把這些頂點按照順時針排個序：
// ⑥ ↓
// 將頂點按照順時針排序，構(gòu)成凸多邊形
std::vector<cv::Point> poly1;
convexHull(outIntersection1, poly1, true, true);
std::vector<cv::Point> poly2;
convexHull(outIntersection2, poly2, true, true);
// ⑥ ↓
看看效果：
最終效果
綠色區(qū)域就是兩張照片覆蓋的部分。
4. 多邊形計算相關(guān)的方法
getPolyIntersections 及相關(guān) callee 方法：
/**
* 計算兩個多邊形的交點。
* 如果沒有交點返回 false。有交點則保存到 outIntersections 中。
*/
bool getPolyIntersections(const std::vector<cv::Point> &poly1, const std::vector<cv::Point> &poly2, std::vector<cv::Point> &outIntersections) {
// 至少是個三角形
if (poly1.size() < 3 || poly2.size() < 3) {
return false;
}
// 計算多邊形交點
long x, y;
for (auto i = 0; i < poly1.size(); i++) {
auto nextIdx1 = (i + 1) % poly1.size();
for (auto j = 0; j < poly2.size(); j++) {
auto nextIdx2 = (j + 1) % poly2.size();
// 計算兩個線段的交點
bool isCross = getSegmentIntersection(
poly1[i],
poly1[nextIdx1],
poly2[j],
poly2[nextIdx2],
x,
y
);
if (isCross) {
outIntersections.emplace_back(x, y);
}
}
}
// 計算多邊形內(nèi)部點
for (const auto &point : poly1) {
if (isPointInPolygon(poly2, point)) {
outIntersections.push_back(point);
}
}
for (const auto &point : poly2) {
if (isPointInPolygon(poly1, point)) {
outIntersections.push_back(point);
}
}
// 如果沒有交點
return !outIntersections.empty();
}
/**
* 獲取兩個線段的交點。
*
* @retun false, 如果沒有交點。
*/
bool getSegmentIntersection(const cv::Point &p1, const cv::Point &p2, const cv::Point &q1, const cv::Point &q2, long &outX, long &outY) {
// 判斷兩條線段是否相交
if (!isSegmentCross(p1.x, p1.y, p2.x, p2.y, q1.x, q1.y, q2.x, q2.y)) {
return false;
}
// 求交點
long left, right;
left = (q2.x - q1.x) * (p1.y - p2.y) - (p2.x - p1.x) * (q1.y - q2.y);
right = (p1.y - q1.y) * (p2.x - p1.x) * (q2.x - q1.x) + q1.x * (q2.y - q1.y) * (p2.x - p1.x) -
p1.x * (p2.y - p1.y) * (q2.x - q1.x);
outX = (int) ((double) right / (double) left);
left = (p1.x - p2.x) * (q2.y - q1.y) - (p2.y - p1.y) * (q1.x - q2.x);
right = p2.y * (p1.x - p2.x) * (q2.y - q1.y) + (q2.x - p2.x) * (q2.y - q1.y) * (p1.y - p2.y) -
q2.y * (q1.x - q2.x) * (p2.y - p1.y);
outY = (int) ((double) right / (double) left);
return true;
}
/**
* 判斷兩條線段是否相交。
*/
bool isSegmentCross(double aX, double aY, double bX, double bY,
double cX, double cY, double dX, double dY) {
// 先做排斥實驗，如果兩條線段組成的兩個矩形沒有重疊，則兩條線段一定沒有相交
if (!isRectCross(aX, aY, bX, bY, cX, cY, dX, dY)) {
return false;
}
// 跨立實驗：
// 已知： 如果 α × β < 0, 則 β 在 α 的順時針方向；同理，>0則在逆時針方向。
// 設(shè)線段l1的端點是 A->B, 線段l2的端點是 C->D
// 先令 α = AC, β = AB 得到 α × β => temp1
// 再令 β = AB, γ = AD 得到 β × γ => temp2
// 如果 temp1 · temp2 > 0， 說明 α -> β -> γ 都是相同方向的，也就是 AB 線段可能穿過了 CD 線段
// 同樣的方式，我們可以算出 CD 線段是否穿過了 AB 線段
// 如果兩個條件都滿足，則 AB線段 與 CD線段 一定相交
// 在計算機中，向量叉積表示為： α × β = α.x · β.y - α.y · β.x
bool temp1 = mayHasIntersection(aX, aY, bX, bY, cX, cY, dX, dY);
if (!temp1) {
return false;
}
return mayHasIntersection(cX, cY, dX, dY, aX, aY, bX, bY);
}
/**
* 從線段1的視角觸發(fā)，判斷線段1是否有可能與線段2相交
*/
bool mayHasIntersection(double aX, double aY, double bX, double bY,
double cX, double cY, double dX, double dY) {
auto acX = cX - aX;
auto acY = cY - aY;
auto abX = bX - aX;
auto abY = bY - aY;
auto adX = dX - aX;
auto adY = dY - aY;
auto temp1 = (acX * abY) - (abX * acY);
auto temp2 = (acX * adY) - (adX * acY);
return temp1 * temp2 > 0;
}
/**
* 判斷兩個矩形(豎邊都平行于Y軸)是否有重疊。
* AB 構(gòu)成一個矩形，CD 構(gòu)成一個矩形。
*/
bool isRectCross(double aX, double aY, double bX, double bY,
double cX, double cY, double dX, double dY) {
// 如果 X 軸沒有相交，肯定沒有重疊
if (std::min(aX, bX) > std::max(cX, dX) /* R2在R1左邊 */ || std::min(cX, dX) > std::max(aX, bX) /* R1在R2左邊 */ ) {
return false;
}
// 如果 Y 軸沒有相交，肯定沒有重疊
if (std::min(aY, bY) > std::max(cY, dY) /* R1在R2上邊 */ || std::min(cY, dY) > std::max(aY, bY) /* R2在R1上邊 */ ) {
return false;
}
return true;
}
/**
* 判斷一個坐標點是否在一個多邊形中。
*/
template<typename Point_t>
static bool isPointInPolygon(const std::vector<Point_t> &poly, const Point_t &pt) {
size_t i, j;
bool c = false;
for (i = 0, j = poly.size() - 1; i < poly.size(); j = i++) {
if ((((poly[i].y <= pt.y) && (pt.y < poly[j].y)) ||
((poly[j].y <= pt.y) && (pt.y < poly[i].y)))
&& (pt.x < (poly[j].x - poly[i].x) * (pt.y - poly[i].y) / (poly[j].y - poly[i].y) + poly[i].x)) {
c = !c;
}
}
return c;
}
希望對你有用~
