6、寻找接缝线
6.1 原理
拼接图像的另一个重要的步骤是找到图像重叠部分内的一条接缝线,该接缝是重叠部分最相似的像素的连线。当确定了接缝线后,在重叠部分,线的一侧只选择该侧的图像部分,线的另一侧只选择这一侧的图像部分,而不是把重叠部分的两幅图像简单融合起来。这么做的目的可以避免图像的模糊及伪像。
目前,常用的寻找接缝线的方法有三种:逐点法、动态规划法和图割法。
逐点法比较简单,它的原理就是重叠部分内的像素距离哪个图像更近,就用哪个图像中的相应像素值。这种方法并不是基于像素值,只是基于距离准则。所以这种方法并不能真正得到最佳的接缝线,寻找的结果比较粗糙原始,仅仅是能够分割重叠部分而已。
动态规划法可以快速的找到最佳的接缝线,并且所需的内存较少,因此该方法很适合在移动终端设备上进行高清晰度的全景拼接。
设图像1和图像2之间有重叠部分,我们需要得到它们之间的最佳接缝,首先定义重叠部分的误差表面函数:
(88)
式中,I1,I2表示两幅图像各自的重叠部分,││,||表示范数。
一般来说,接缝线有三个限制条件:一是如果重叠区域是宽大于高,则接缝是横向走向的,而如果重叠区域是宽小于高,则接缝是纵向走向的,也就是要保证接缝线要有一定的长度;二是如果是横向的接缝,则不允许有绝对垂直的接缝线,而如果是纵向的接缝,则不允许有绝对水平的接缝线;三是重叠区域是矩形,接缝线从矩形的一边出发,必须到达与该边平行的另一边结束。
假设重叠区域的宽小于高,则接缝线是垂直的,我们需要横向遍历e值,并且计算所有可能到达当前像素(i, j)的路径的累积最小误差E:
(89)
式中,min中的三项内容分别表示当前像素与其左上侧、上侧和右上侧的E。在E中最后一行的最小值表明已经到达了最小垂直路径的尽头,那么我们就可以追溯,从而得到最佳路径,即接缝线。
同理,重叠区域的宽大于高也能够做出类似的计算过程。
我们可以看出,动态规划算法非常适用于对式89的求解,这也是该寻找最佳接缝线方法名称的由来。
对于彩色图像,计算式88可以有两种方法:直接法和梯度法。设p和q是两个相邻像素,直接法的计算公式为:
(90)
梯度法的计算公式为:
(91)
式中,G表示灰度图像的p到q方向上的梯度。
最后一种方法是图割法。它的主要思想是将图像看成是一个有向图G={ V, E},V表示节点,E表示边。节点分为两种:普通节点和终端节点。图像的像素可以作为普通节点,图像的聚类结果可以看成是终端节点,即某一聚类中的普通节点会汇聚于一个终端节点,一般图只有两个终端节点:s和t,也就是说一般图割法只能够把图像分割成两类S和T。边指的是节点之间的连接,对于普通节点,只有相邻节点可以连接,该边被称为n-links,而所有普通节点都可与属于它们的终端节点之间有连接,该边被称为t-links。用最大流图割法可以得到哪些普通节点属于终端节点s,那么剩下的普通节点就会属于t。
Opencv中所实现的最大流图割算法是基于Boykov和Kolmogorov所提出的算法。该算法是以s和t为根,以得到的两棵无不相交的搜索树S和T为目标,它迭代的重复以下三个阶段:
生长阶段:搜索树S和T的生长,直到两棵树生长到了一起,形成了从s到t的通路(即为流);
增广阶段:该通路被增广,使相连的两棵树又重新分开;
收养阶段:搜索树S和T收养孤立的节点(称为孤儿)。
在生长阶段,如果找不到任何一条从s到t的通路,则迭代结束。在增广阶段,开始孤儿集合是空集,而该阶段结束后,可能会出现一些孤儿。在该阶段,我们需要得到通路上最小的容量值(也称为边的权值),然后所有的容量值都减去该值(即通过流量更新容量值),则该通路上至少会有一个权值为0的边(称为饱和边),假设饱和边上的两个节点为p和q,并且在正向通路中是由p流向q,如果p和q都属于树S,则q为孤儿,如果p和q都属于树T,则p为孤儿。在收养阶段,检查孤儿集合中的所有孤儿的邻域像素,连接那些它们的边的容量值不为0的像素,如果最终可以从该孤儿连通到s或t,则该孤儿就属于树S或树T,否则由该孤儿组成的树就属于自由节点集,在下次迭代时的生长阶段,树S和树T就通过连接自由节点集中的节点来使两棵树生长发育的。
如何用最大流图割法寻找接缝线呢?接缝线把重叠区域分割成了两个区域,那么我们通过最大流图割法得到树S和树T,即得到了被分割的两个区域,则区域之间的接缝线自然就可以得到。
普通节点之间的初始容量值(即边的权值)也可以通过直接法和梯度法计算得到,它们的计算公式分别与式90和式91相同。终端节点与普通节点之间的初始容量值可以定义为一个很大的值,当重叠区域中的节点只属于图像1或图像2时,节点与相应的终端节点的边赋予该值,其他情况都置为0,表示它们与终端节点无直接连接。
6.2 源码
SeamFinder类是寻找缝的基类,该类内有一个最重要的虚函数是find,它的作用是寻找缝:
virtual void find(const std::vector<Mat> &src, const std::vector<Point> &corners,
std::vector<Mat> &masks) = 0;
它的三个输入参数分别为:src表示待拼接的图像,corners表示待拼接图像在最终的全景图像坐标系内的左上角的坐标位置,masks输入时表示上一步得到的映射变换掩码,输出时表示接缝掩码,也就是该参数是要更新输出的,输出时如果像素值更新为0了,则表示该像素是在接缝的另一侧,拼接时是不需要该像素的。
基于SeamFinder类有四个子类,分别对应于四种不同的寻找缝的算法——NoSeamFinder(无需寻找缝)、PairwiseSeamFinder(逐点寻找缝算法)、DpSeamFinder(动态规划寻找缝算法)和GraphCutSeamFinder(图割寻找缝算法)。
我们先来介绍PairwiseSeamFinder类。逐点寻找缝算法又可分为不同的方法,Opencv只实现了Voronoi图方法,实现该方法的类为VoronoiSeamFinder类,它是PairwiseSeamFinder类的子类。VoronoiSeamFinder类的find函数为:
void VoronoiSeamFinder::find(const vector<Size> &sizes, const vector<Point> &corners,
vector<Mat> &masks)
{
LOGLN("Finding seams...");
if (sizes.size() == 0) //没有图像,则返回
return;
#if ENABLE_LOG
int64 t = getTickCount(); //计时
#endif
//参数赋值
sizes_ = sizes;
corners_ = corners;
masks_ = masks;
run(); //调用run函数
//计时输出
LOGLN("Finding seams, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
}
VoronoiSeamFinder类没有实现run函数,因此这里的find函数调用的是它的父类PairwiseSeamFinder类的run函数:
void PairwiseSeamFinder::run()
{
//遍历图像对
for (size_t i = 0; i < sizes_.size() - 1; ++i)
{
for (size_t j = i + 1; j < sizes_.size(); ++j)
{
Rect roi; //表示图像间的重叠部分
//如果图像i和图j之间有重叠的区域,则执行findInPair函数
if (overlapRoi(corners_[i], corners_[j], sizes_[i], sizes_[j], roi))
findInPair(i, j, roi);
}
}
}
findInPair函数表示寻找两幅图像之间的接缝:
void VoronoiSeamFinder::findInPair(size_t first, size_t second, Rect roi)
//first和second表示两幅图像的尺寸大小,roi表示图像间重叠部分
{
const int gap = 10;
//分别表示图像first和图像second的重叠部分的掩码,它们要比roi四周多gap个距离
Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
Size img1 = sizes_[first], img2 = sizes_[second]; //得到两幅图像的尺寸
Mat mask1 = masks_[first], mask2 = masks_[second]; //得到两幅图像的掩码
Point tl1 = corners_[first], tl2 = corners_[second]; //得到两幅图像的左上角坐标
// Cut submasks with some gap
//遍历重叠部分,得到submask1和submask2
for (int y = -gap; y < roi.height + gap; ++y)
{
for (int x = -gap; x < roi.width + gap; ++x)
{
//y1和x1表示重叠部分相对于图像first的坐标
int y1 = roi.y - tl1.y + y;
int x1 = roi.x - tl1.x + x;
//把图像first的掩码赋值给submask1
//表示当前坐标在图像first内
if (y1 >= 0 && x1 >= 0 && y1 < img1.height && x1 < img1.width)
submask1.at<uchar>(y + gap, x + gap) = mask1.at<uchar>(y1, x1);
else //表示当前坐标不在图像first内
submask1.at<uchar>(y + gap, x + gap) = 0;
//y2和x2表示重叠部分相对于图像second的坐标
int y2 = roi.y - tl2.y + y;
int x2 = roi.x - tl2.x + x;
//把图像second的掩码赋值给submask2
//表示当前坐标在图像second内
if (y2 >= 0 && x2 >= 0 && y2 < img2.height && x2 < img2.width)
submask2.at<uchar>(y + gap, x + gap) = mask2.at<uchar>(y2, x2);
else //表示当前坐标不在图像second内
submask2.at<uchar>(y + gap, x + gap) = 0;
}
}
//collision表示submask1和submask2交集,roi区域可能有些区域是无效的部分(由于映射的原因),而collision(像素值为1的部分)则是图像first和图像second的真正有意义部分的重叠区域,我们是需要在collision区域内寻找接缝的
Mat collision = (submask1 != 0) & (submask2 != 0);
//复制submask1和submask2为unique1和unique2,并且把collision部分的像素清零
Mat unique1 = submask1.clone(); unique1.setTo(0, collision);
Mat unique2 = submask2.clone(); unique2.setTo(0, collision);
Mat dist1, dist2; //表示距离矩阵
//分别得到unique1和unique2内为0的像素值与最近的非0值之间的L1距离,显然为0的区域肯定要比collision大
distanceTransform(unique1 == 0, dist1, CV_DIST_L1, 3);
distanceTransform(unique2 == 0, dist2, CV_DIST_L1, 3);
//得到缝矩阵,dist1 < dist2成立,说明离first图像近,则seam内相应的像素为1,否则为0
Mat seam = dist1 < dist2;
//遍历roi内的所有像素
for (int y = 0; y < roi.height; ++y)
{
for (int x = 0; x < roi.width; ++x)
{
//如果seam内的像素值为1,说明collision区域内的相应像素离图像first近,所以该点的值选择图像first的像素值,因此需把mask2清零,反之把mask1清零。也就是在重叠部分,哪些区域应该属于图像first,哪些区域属于图像second。
if (seam.at<uchar>(y + gap, x + gap))
mask2.at<uchar>(roi.y - tl2.y + y, roi.x - tl2.x + x) = 0;
else
mask1.at<uchar>(roi.y - tl1.y + y, roi.x - tl1.x + x) = 0;
}
}
}
DpSeamFinder类利用动态规划算法寻找缝,它的构造函数为:
DpSeamFinder::DpSeamFinder(CostFunction costFunc) : costFunc_(costFunc) {}
其中CostFunction为:
enum CostFunction { COLOR, COLOR_GRAD };
它表示动态规划寻找接缝的两种不同方法,直接法和梯度法,默认是直接法。
DpSeamFinder类的find函数为:
void DpSeamFinder::find(const vector<Mat> &src, const vector<Point> &corners, vector<Mat> &masks)
{
LOGLN("Finding seams...");
#if ENABLE_LOG
int64 t = getTickCount(); //用于计时
#endif
if (src.size() == 0) //没有图像,则返回
return;
vector<pair<size_t, size_t> > pairs; //表示图像对
for (size_t i = 0; i+1 < src.size(); ++i)
for (size_t j = i+1; j < src.size(); ++j)
pairs.push_back(make_pair(i, j)); //每两幅图像组成一个图像对
//根据图像位置(基于左上角坐标),对图像对进行排序
sort(pairs.begin(), pairs.end(), ImagePairLess(src, corners));
reverse(pairs.begin(), pairs.end()); //逆序
for (size_t i = 0; i < pairs.size(); ++i) //遍历图像对,
{
size_t i0 = pairs[i].first, i1 = pairs[i].second;
//寻找i0和i1这两幅图像的接缝,该函数在后面有介绍
process(src[i0], src[i1], corners[i0], corners[i1], masks[i0], masks[i1]);
}
LOGLN("Finding seams, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
}
寻找两幅图像间的接缝:
void DpSeamFinder::process(
const Mat &image1, const Mat &image2, Point tl1, Point tl2,
Mat &mask1, Mat &mask2)
//image1和image2分别表示待处理的两幅图像
//tl1和tl2分别表示这两幅图像的左上角坐标
//mask1和mask2分别表示这两幅图像的掩码
{
CV_Assert(image1.size() == mask1.size()); //确保图像与掩码尺寸相同
CV_Assert(image2.size() == mask2.size());
//intersectTl和intersectBr分别表示image1和image2的交集区域(即重叠部分)的左上角和右下角坐标
Point intersectTl(std::max(tl1.x, tl2.x), std::max(tl1.y, tl2.y));
Point intersectBr(std::min(tl1.x + image1.cols, tl2.x + image2.cols),
std::min(tl1.y + image1.rows, tl2.y + image2.rows));
//如果if条件成立,则说明image1和image2没有重叠区域,因此退出该函数
if (intersectTl.x >= intersectBr.x || intersectTl.y >= intersectBr.y)
return; // there are no conflicts
//unionTl_和unionBr_分别表示image1和image2的并集区域(即外围最大部分的矩形)的左上角和右下角坐标
unionTl_ = Point(std::min(tl1.x, tl2.x), std::min(tl1.y, tl2.y));
unionBr_ = Point(std::max(tl1.x + image1.cols, tl2.x + image2.cols),
std::max(tl1.y + image1.rows, tl2.y + image2.rows));
//得到image1和image2的并集区域的尺寸大小
unionSize_ = Size(unionBr_.x - unionTl_.x, unionBr_.y - unionTl_.y);
//定义两幅图像的掩码大小,先清零
mask1_ = Mat::zeros(unionSize_, CV_8U);
mask2_ = Mat::zeros(unionSize_, CV_8U);
//tmp为mask1_中image1的区域
Mat tmp = mask1_(Rect(tl1.x - unionTl_.x, tl1.y - unionTl_.y, mask1.cols, mask1.rows));
mask1.copyTo(tmp); //把image1的掩码mask1赋值给mask1_
//tmp为mask2_中image2的区域
tmp = mask2_(Rect(tl2.x - unionTl_.x, tl2.y - unionTl_.y, mask2.cols, mask2.rows));
mask2.copyTo(tmp); //把image2的掩码mask2赋值给mask2_
// find both images contour masks
//在并集中分别得到边界掩码contour1mask_和contour2mask_
contour1mask_ = Mat::zeros(unionSize_, CV_8U); //先清零
contour2mask_ = Mat::zeros(unionSize_, CV_8U);
//遍历并集区域
for (int y = 0; y < unionSize_.height; ++y)
{
for (int x = 0; x < unionSize_.width; ++x)
{
if (mask1_(y, x) &&
((x == 0 || !mask1_(y, x-1)) || (x == unionSize_.width-1 || !mask1_(y, x+1)) ||
(y == 0 || !mask1_(y-1, x)) || (y == unionSize_.height-1 || !mask1_(y+1, x))))
{
contour1mask_(y, x) = 255; //得到image1的边界掩码
}
if (mask2_(y, x) &&
(