图片重叠区域识别有很多应用场景，例如全景照片的合成等。这篇文章介绍一种图片重叠区域识别的方法。
1. 整体流程
识别重叠区的步骤包含：
- 查找两张图片的特征点；
- 匹配特征点；
- 计算单应矩阵；
- 基于单应矩阵，计算第二张图片的四个顶点在图片1中的位置；
- 计算 (4) 中得到的位置与图片边缘的交点，得到重叠区域的顶点；
- 同理得到图片1的顶点在图片2中覆盖区域的顶点。
来看代码和效果：
2. 查找和匹配特征点
找特征点的方式有很多，这里用的是 SIFT，也可以用 SURF、ORB 都行。目的是为了计算出单应矩阵。
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);
// ① ↓
经过这两步我们就能得到匹配的特征点了，但这些匹配对有很多不正确的：
匹配后的特征点
通过汉明距离过滤掉部分匹配对：
// ① ↓
// 获取匹配对汉明距离的最小值和最大值
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);
}
}
// ② ↓
过滤后的匹配对：
挑选后的匹配对
不同的照片这一步得到的匹配对数量不一样，通常变化不大的场景留下的匹配对会更多。
3. 单应矩阵和重叠区域的计算
单应矩阵约束了同一3D空间点在两个像素平面的2D齐次坐标。有了单应矩阵我们就能计算图片2的顶点在图片1中应该在什么位置。
// ② ↓
// 找出关键点对应的像素坐标
std::vector<cv::Point2f> points1, points2;
for (auto &match : goodMatches) {
points1.push_back(keyPoint2[match.trainIdx].pt);
points2.push_back(keyPoint1[match.queryIdx].pt);
}
// 计算单应矩阵
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 的四个顶点，根据单应矩阵，计算出这四个顶点在 图片1 中的位置
for (auto& point : pointsB) {
cv::Point temp;
// 这个方法在下面👇
transformPoint(H, point, temp);
vertices2AtPhoto1.push_back(temp);
}
// ④ ↓
其中 transformPoint 为：
/**
* 计算坐标 in 经过 单应矩阵 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中的位置，但是这个位置可能已经超出图片1的视野外了，所以想要计算重叠区域，还需要计算图2变换后的顶点与图一的交点。交点围成的多边形就是重叠区域。这是个纯数学问题：
// ④ ↓
// 计算 图像B 的顶点在图像A中 和 图像A边框 的交点
std::vector<cv::Point> outIntersection1;
getPolyIntersections(pointsB, vertices2AtPhoto1, outIntersection1);
// ⑤ ↓
其中 getPolyIntersections 及相关 callee 方法我放在了文章最后。
同样，我们计算出图1在图2中的多边形顶点：
// ⑤ ↓
// 用同样的方式，计算出图片1的四个顶点，在图片2中与图片2的交点。
auto hInverses = H.inv();
std::vector<cv::Point> outIntersection2;
// 遍历 图片2 的四个顶点，根据单应矩阵，计算出这四个顶点在 图片1 中的位置
for (auto& point : vertex) {
cv::Point temp;
transformPoint(hInverses, point, temp);
outIntersection2.push_back(temp);
}
// ⑥ ↓
为了更好的显示和后续覆盖度的计算等，我们可以把这些顶点按照顺时针排个序：
// ⑥ ↓
// 将顶点按照顺时针排序，构成凸多边形
std::vector<cv::Point> poly1;
convexHull(outIntersection1, poly1, true, true);
std::vector<cv::Point> poly2;
convexHull(outIntersection2, poly2, true, true);
// ⑥ ↓
看看效果：
最终效果
绿色区域就是两张照片覆盖的部分。
4. 多边形计算相关的方法
getPolyIntersections 及相关 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);
}
}
}
// 计算多边形内部点
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则在逆时针方向。
// 设线段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的视角触发，判断线段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 构成一个矩形，CD 构成一个矩形。
*/
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;
}
希望对你有用~
