全景图像拼接
- 全景图像拼接的手动实现
- 环境： python3.6 + opencv3.4.2.16
## 示例图片
本次实验使用的图像拼接素材为以下三张图像：
https://andreame.com/2019/11/12/stitch.html
本次实验的目标为，将此三张图像进行圆柱面投影并进行全景拼接
opencv内置实现
首先，opencv已经内置了stitch这一个类，包装好了全景图像拼接的所能用到的方法。
基础使用相当简单，核心代码两行实现
import cv2
from cv2 import Stitcher
def stitcher(images)
# images是待拼接图像的list
Stitcher = cv2.Stitcher.create(cv2.Stitcher_PANORAMA)
(status, pano) = Stitcher.stitch(images)
return panoPython
Copy
并且默认的拼接效果还相当不错
所以仅仅是实现完成图像拼接功能的话，以上两行代码就已经可以解决了。
分析
当然，仅仅学会调用是远远不够的，最好还是自己尝试手动实现，很可惜opencv是调用的底层cpp库，而非python实现，直接修改源码的想法并不能实施，封装好的stitch类的Pipeline如下
opencv文档指出，这个类主要包含以下几个模块：
- Features Finding and Images Matching 功能查找和图像匹配
- Rotation Estimation 旋转估计
- Autocalibration 自动校准
- Images Warping 图像变形
- Seam Estimation 接缝估计
- Exposure Compensation 曝光补偿
- Image Blenders 图像拼接
着一些拼接过程速度较慢，对于普通的全景图像拼接，并不需要如此繁琐的步骤，因此我们可以自行设计模块组成。
设计与实现
代码实现主要是参考了这一篇博客与这一篇教程 主要步骤如下：
- 图像的圆柱投影变换
- 计算待拼接图像的特征点以及描述点
- 计算图像之间的特征描述位置距离
- 筛选最优的特征点
- 使用RANSAC计算单应矩阵
- Warping
- 接缝融合
- Stitching
- 图像裁剪
圆柱投影
简单地说，将原图像上某一个像素(x,y)投影到一个焦距为f的圆柱面上则
\begin{array}{l}{\text { 1. } \theta=\arctan \left(\frac{x}{f}\right)} \\ {\text { 2. } x^{\prime}=f \times \theta} \\ {\text { 3. } \frac{y}{\frac{f}{\cos (\theta)}}=\frac{y^{\prime}}{f}}\end{array} 1. θ=arctan(fx) 2. x′=f×θ 3. cos(θ)fy=fy′
通过上式我们可以做到$(x,y)$到$(x^{\prime},y^{\prime})$的变化，代码如下：
def cylindricalProjection(img) :
rows = img.shape[0]
cols = img.shape[1]
#f = cols / (2 * math.tan(np.pi / 8))
result = np.zeros_like(img)
center_x = int(cols / 2)
center_y = int(rows / 2)
alpha = math.pi / 4
f = cols / (2 * math.tan(alpha / 2))
for y in range(rows):
for x in range(cols):
theta = math.atan((x- center_x )/ f)
point_x = int(f * math.tan( (x-center_x) / f) + center_x)
point_y = int( (y-center_y) / math.cos(theta) + center_y)
if point_x >= cols or point_x < 0 or point_y >= rows or point_y < 0:
pass
else:
result[y , x, :] = img[point_y , point_x ,:]
return resultPython
Copy
变换结果如下
右图是经过左图经过圆柱面投影之后的结果。
当然，为了保证黑边不影响图像拼接的效果，我们还需要提前把黑边去掉。如何去除黑边，我们会在之后提到。
特征点提取
特征点提取
特征点提取可以使用opencv内置的函数detectAndCompute()实现。
需要注明的是，在opencv-python4以上的版本中，由于版权问题，此函数已不再开放，需要将版本退回到opencv3。笔者使用的环境为
- python3.6
- opencv-contrib-python 3.4.2.16
- opencv-python 3.4.2.16
函数调用如下：
def getSURFFeatures(im)
gray = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
surf = cv2.xfeatures2d.SURF_create()
kp, des = surf.detectAndCompute(gray, None)
#sift = cv2.xfeatures2d.SIFT_create()
#kp, des = sift.detectAndCompute(im, None)
cv2.imwrite('imageWithKeypoints.png',cv2.drawKeypoints(im,kp,None))
return {'kp':kp, 'des':des}Python
Copy
这里可以选用sift方法或者opencv内置的surf方法。下图展示了这些图3里的特征点
特征点距离计算
我们得到了两张图片的特征点之后，我们下一步要做的就是找到这两张图片之间特征点的对应关系。对于将两张图片拼接成一张大图像，我们需要找到他们之间的重合点。利用这些个重合点可以确定如何将第一张图像进行怎样的旋转缩放变形，使得能与下一张图像吻合拼接起来。这就是找到一种对应关系。
这里可以使用opencv提供的两种方法去匹配图像FLANN或BFMatcher方法。
下面是具体的实现代码：
# matches = self.flann.knnMatch(
# imageSet2['des'],
# imageSet1['des'],
# k=2
# )
match = cv2.BFMatcher()
matches = match.knnMatch(
imageSet2['des'],
imageSet1['des'],
k=2
)Python
Copy
其中k=2，这个参数是指让matches给出两个最佳匹配。
draw_params = dict(matchColor = (0,255,0), # draw matches in green color
singlePointColor = None,
flags = 2)
img3 = cv2.drawMatchesKnn(i1,imageSet1['kp'],i2,imageSet2['kp'],matches,None,**draw_params)
cv2.imwrite("correspondences.png", img3)
print(len(matches))Python
Copy
绘制下效果，下图展示了两张图像之间的特征点匹配。
显然，特征点数量过多，特征点数量为615
筛选最优匹配特征点
通常在图像中，特征点存在于图像中的多个位置，因此需要进行进一步筛选。
good = []
for i, (m, n) in enumerate(matches):
if m.distance < 0.3*n.distance:
good.append((m.trainIdx, m.queryIdx))Python
Copy
这样的话，我们就得到了两张图像之间的最佳匹配，下一步则是计算单应矩阵。
当这个ratio设为0.3时，筛选的效果如下图所示
此时特征点数量为39
单应矩阵的计算
单应矩阵的作用是依据最佳匹配点，估计两幅图像中的相对方向变换。这里可以使用cv内置的RANSAC方法进行求解：
H, s = cv2.findHomography(matchedPointsCurrent, matchedPointsPrev, cv2.RANSAC, 4)Python
Copy
得到的$H$则是要求解的单应矩阵。
I_x = H \times I_yIx=H×Iy
得到的单应性矩阵如下式
\begin{bmatrix} h_\text{11} & h_\text{12} & h_\text{13} \\ h_\text{21} & h_\text{22} & h_\text{23} \\ h_\text{31} & h_\text{32} & h_\text{33} \end{bmatrix}⎣⎡h11h21h31h12h22h32h13h23h33⎦⎤
Wrap
在上一步中，我们得到了一个单应矩阵，下一步，则是进行图象的缝合。
首先，我们得到了两者之间的单应矩阵，我们就可以了解到，从第一张图象的角度去看第二张图像是什么样的了。我们就需要将图像进行变化，转换到一个新空间。那么接下来就需要进行一个翘曲变换的过程。
翘曲变换包括
- 平面变换
- 圆柱面变换
- 球形变换
我们可以简单调用以下函数来完成。
warped_image = cv2.warpPerspective(image, homography_matrix, dimension_of_warped_image)Python
Copy
以上函数中，image代表了待变换图像，homography_matrix就是我们求得的单应矩阵，而dimension_of_warped_image是变换之后，图像的尺寸。
因此，为了做出翘曲变换，我们还需要求出翘曲变换之后的尺寸，为拼接做准备。
H = self.matcher_obj.match(a, b, 'left')
xh = np.linalg.inv(H)
ds = np.dot(xh, np.array([a.shape[1], a.shape[0], 1]))
ds = ds / ds[-1] #
f1 = np.dot(xh, np.array([0, 0, 1]))
f1 = f1 / f1[-1]
xh[0][-1] += abs(f1[0])
xh[1][-1] += abs(f1[1]) #
ds = np.dot(xh, np.array([a.shape[1], a.shape[0], 1]))
offsety = abs(int(f1[1]))
offsetx = abs(int(f1[0]))
dsize = (int(ds[0]) + offsetx, int(ds[1]) + offsety)Python
Copy
最后得出的可以求得翘曲变换之后的尺寸。
我们有一个单应矩阵$H$ 。如果说每张图片的起始坐标为$(0,0)$ 结束的坐标为$(r_e,c_e)$ ，那么我们可以通过，这样的变化得到翘曲变换后的图像尺寸：
起始点计算$Point_{strat} := H \times (0,0)$ 一直到终点计算 $point_{end} := H \times (r-e,c_e)$。
当然，下一步的拼接可以简单地用tmp[offsety:b.shape[0]+offsety, offsetx:b.shape[1]+offsetx] = b来完成，但这样会使得接缝非常明显，因此我们还需要进行调整。
效果如下：
可以看到非常明显的接缝，显然这么简单地处理并不能达到理想的效果。
接缝融合
那么如何消除接缝呢，这里采用加权融合的方法。
思路是计算出两张图像重叠的部分，对于有重叠的部分，采用像素逐渐过渡的方式去消弭这种肉眼上能观测出的突兀变化。
那么可以考虑使用加权融合的方式去解决，即让拼接部分的像素点为两张图像同坐标的像素点的加权和：
point_{result} = \alpha \times point_{a} + (1 - \alpha) \times point_{b}pointresult=α×pointa+(1−α)×pointb
公式中的$\alpha$如何确定，这里是用距离来计算这个权重，图像从重合边界到图像边界处，这个权重平滑地从1将为0。
def seamEstimation(self,tmp,leftimage,rightimage,offsetx,offsety):
alpha = 0.0
processwidth = leftimage.shape[1] - offsetx
print("initial_width",processwidth)
for x in range(offsetx,tmp.shape[1]):
test_y = int((offsety + leftimage.shape[0])/2)
if np.array_equal(tmp[test_y,x],np.array([0,0,0])):
#print(tmp[test_y,x])
processwidth = x - offsetx
break
print("now_width",processwidth)
x_max = np.minimum(offsetx+rightimage.shape[1],tmp.shape[1])
y_max = np.minimum(offsety+rightimage.shape[0],tmp.shape[0])
for x in range(offsetx,x_max):
for y in range(offsety,y_max):
rx = x - offsetx
ry = y - offsety
if not np.array_equal(tmp[y,x],np.array([0,0,0])) and not np.array_equal(rightimage[ry,rx],np.array([0,0,0])):
alpha = np.minimum(rx/processwidth,1.0)
#alpha = rx/processwidth
tmp[y,x] = alpha*rightimage[ry,rx]+(1-alpha)*tmp[y,x]
elif np.array_equal(tmp[y,x],np.array([0,0,0])):
tmp[y,x] = rightimage[ry,rx]
elif np.array_equal(rightimage[ry,rx],np.array([0,0,0])):
tmp[y,x] = tmp[y,x]
else:
print("Some problems happen!")
return tmpPython
Copy
然后多张图像连续进行拼接，可以得到拼接好的图像:
整体效果还可以，但是明显拼接部位出现了重影，当然这是由于像素融合所造成的必然结果，这一点需要在之后进行改进。
MultipleImagesStitching
通过上述描述，我们已经完成了图像拼接的基本工作，那么接下来仅仅需要不断重复上述的俩俩拼接工作即可。
需要注意的是，如果采取平面变化的话，需要以中间的一张图片为基线，进行左右拼接，以免偏差过大。并且当图片数量过大时，则需要进行分批拼接。
但这里采用的是圆柱面拼接，依次拼接不会造成太过严重的后果。拼接的代码如下：
def shift(self):
a = self.images[0]
for b in self.images[1:]:
H = self.matcher_obj.match(a, b, 'left') # 特征点匹配
xh = np.linalg.inv(H)
ds = np.dot(xh, np.array([a.shape[1], a.shape[0], 1]))
ds = ds / ds[-1] #
f1 = np.dot(xh, np.array([0, 0, 1]))
f1 = f1 / f1[-1]
xh[0][-1] += abs(f1[0])
xh[1][-1] += abs(f1[1]) #
ds = np.dot(xh, np.array([a.shape[1], a.shape[0], 1]))
offsety = abs(int(f1[1])) # y偏移量 需要了解单应性矩阵的作用
offsetx = abs(int(f1[0])) # x偏移量
dsize = (int(ds[0]) + offsetx, int(ds[1]) + offsety) # 图片大小统计
tmp = cv2.warpPerspective(a, xh, dsize)
tmp = self.seamEstimation(tmp,a,b,offsetx,offsety)
a = tmp # 为循环做准备
return aPython
Copy
图像裁剪
当然，我们此时得到的图像还有大量的黑框，我们需要通过图像操作，进行裁剪。
将图像转化为灰度图
img = cv2.medianBlur(image,3) #中值滤波
b=cv2.threshold(img,0,255,cv2.THRESH_BINARY)
binary_image=b[1] #二值图--具有三通道
binary_image=cv2.cvtColor(binary_image,cv2.COLOR_BGR2GRAY)Python
Copy
这里简单运用图像操作将图像转化为二值图。
除此之外我们还需要继续外展10个像素，确保他的边界全部露出，为下一步服务。
计算最大轮廓
ret, thresh = cv2.threshold(binary_image, 0, 255, cv2.THRESH_BINARY)
binary ,cnts, hierarchy = cv2.findContours(thresh.copy(), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
cnt = max(cnts, key=cv2.contourArea)Python
Copy
使用opencv的内置函数来计算这个图像的最大轮廓。
绘制最大外接矩形框
mask = np.zeros(thresh.shape, dtype="uint8")
x, y, w, h = cv2.boundingRect(cnt)Python
Copy
连续腐蚀操作
我们将初始边界框定在最大外接矩形框上，然后以此为起点，不断地进行服饰操作，缩减所选框的面积，知道，我们的所选框内不在含有背景像素位置，这样我们就找到了我们所需的边界框。$(x,y,w,h)$
下图则是裁剪后的拼接图像。
以下为去除边框的代码：
def ExtractImage(image):
# 去除边框，提取图像内容
plt.figure(num='ExtractImage')
image = cv2.copyMakeBorder(image, 10, 10, 10, 10, cv2.BORDER_CONSTANT, (0, 0, 0))
binary_image = getBinaryImage(image)
cv2.imwrite("binary_image.png", binary_image)
ret, thresh = cv2.threshold(binary_image, 0, 255, cv2.THRESH_BINARY)
binary ,cnts, hierarchy = cv2.findContours(thresh.copy(), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
cnt = max(cnts, key=cv2.contourArea) # 获取最大轮廓
mask = np.zeros(thresh.shape, dtype="uint8")
x, y, w, h = cv2.boundingRect(cnt)
# 绘制最大外接矩形框（内部填充）
cv2.rectangle(mask, (x, y), (x + w, y + h), 255, -1)
cv2.imwrite("mask.png", mask)
minRect = mask.copy()
sub = mask.copy()
print(sub.shape[0] * sub.shape[1])
# 连续腐蚀操作，直到sub中不再有前景像素
while cv2.countNonZero(sub) > 0:
minRect = cv2.erode(minRect, None)
sub = cv2.subtract(minRect, thresh)
binary ,cnts, hierarchy = cv2.findContours(minRect.copy(), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
cnt = max(cnts, key=cv2.contourArea)
x, y, w, h = cv2.boundingRect(cnt)
image = image[y:y + h, x:x + w]
return imagePython
Copy
孔洞填充
上一步中，我们设定腐蚀操作的循环条件是“候选框内不得有背景像素”，那么如何判断背景像素呢？可以简单地认为像素值为0即为背景像素。但如果原图像中，或者其他操作中会导致图像内容中出现像素值为0的像素，如：
这两张图像拼接之后会得到的图转化为的二值图，内有空洞，贸然裁剪会减去上半部分，
这样的图像，直接进行裁剪会导致以下结果：
因此还需要先进性孔洞填充，确保图像内没有孔洞，才可进行下一步的裁剪
填充过后，裁剪结果：
最终结果
拼接前
拼接后
缩放
Time cost :5.6752543449401855
