【计算机视觉】图像映射和全景拼接
1. 图像映射
图像映射流程
1.提取特征点，生成描述符
2.特征匹配
3.根据图像变换特点，选取合适的变换类型
4.根据单应性变换等方法计算变换结构
5.采用正向或逆向映射，利用插值方式实现图像映射变换
图像变换类型
1.刚体变换：只改变物体位置，不改变物体形状（如平移、旋转）
2.仿射变换：改变物体位置和形状，但是保持“平直性”
3.投影变换：彻底改变物体位置和形状
1.1 单应性变换
单应变换是将一个平面中的点映射到另一个平面的二维投影变换。实现单应性变换的关键是得到单应性矩阵，对于不同的变换类型，需要对应数量的对应点对才能得到对应的单应性矩阵（求解矩阵中的未知数）。
不同变换需要相应的对应点对（一个点对分别得到x和y两个方程）
平移变换：2个自由度（需要一个点对）
相似变换：4个自由度（需要两个点对）
仿射变换：6个自由度（需要三个点对）
投影变换：8个自由度（需要四个点对）
求解出单应性矩阵后，我们可以将其带入图像的所有像素中，得到单应性变换后的图像。求解单应矩阵的过程如下：
（1） 其中x,y为原图像的像素坐标，x’,y’为变换后图像的像素坐标
（2） 令w=1，目的：进行点的归一化
（3） 将其转化为Ah = 0的形式
（4） 当我们凑齐Ah = 0的公式后，A是一个对应点对二倍数量的行数的一个矩阵（比如要实现投影变换需要4个点对，这个A矩阵就有2×4=8行）。所以我们根据不同的图像变换类型设置不同的点对数n，A矩阵的维度为2n x 9，h矩阵的维度为9，0矩阵的维度为2n，最后根据最小二乘法进行求解。
（5） 代码验证
import numpy as np
from pylab import *
from PIL import Image
def H_from_points(fp,tp):
""" Find homography H, such that fp is mapped to tp
using the linear DLT method. Points are conditioned
automatically. """
if fp.shape != tp.shape:
raise RuntimeError('number of points do not match')
# condition points (important for numerical reasons)
# --from points--
m = mean(fp[:2], axis=1)
maxstd = max(std(fp[:2], axis=1)) + 1e-9
C1 = diag([1/maxstd, 1/maxstd, 1])
C1[0][2] = -m[0]/maxstd
C1[1][2] = -m[1]/maxstd
fp = dot(C1,fp)
# --to points--
m = mean(tp[:2], axis=1)
maxstd = max(std(tp[:2], axis=1)) + 1e-9
C2 = diag([1/maxstd, 1/maxstd, 1])
C2[0][2] = -m[0]/maxstd
C2[1][2] = -m[1]/maxstd
tp = dot(C2,tp)
# create matrix for linear method, 2 rows for each correspondence pair
nbr_correspondences = fp.shape[1]
A = zeros((2*nbr_correspondences,9))
for i in range(nbr_correspondences):
A[2*i] = [-fp[0][i],-fp[1][i],-1,0,0,0,
tp[0][i]*fp[0][i],tp[0][i]*fp[1][i],tp[0][i]]
A[2*i+1] = [0,0,0,-fp[0][i],-fp[1][i],-1,
tp[1][i]*fp[0][i],tp[1][i]*fp[1][i],tp[1][i]]
U,S,V = linalg.svd(A)
H = V[8].reshape((3,3))
# decondition
H = dot(linalg.inv(C2),dot(H,C1))
# normalize and return
return H / H[2,2]
if __name__ == "__main__":
fp = np.random.randint(0, 500, size=(3, 3))
tp = np.random.randint(0, 500, size=(3, 3))
result = H_from_points(fp,tp)
print(result)
（6） 结果验证，根据代码求解出的单应性矩阵满足八个自由度的求解
[[-3.69744393e+00 5.73332427e+00 -3.60193467e+02]
[-1.31628505e+00 2.40124692e+00 -1.83652536e+03]
[-2.12511707e-03 3.04655705e-03 1.00000000e+00]]
（7） 可视化
from scipy import ndimage
from PIL import Image
from pylab import *
im = array(Image.open('sun.jpg').convert('L'))
H = array([[1.4,0.05,-100],[0.05,1.5,-100],[0,0,1]])
im2 = ndimage.affine_transform(im,H[:2,:2],(H[0,2],H[1,2]))
figure()
gray()
subplot(121)
axis('off')
imshow(im)
subplot(122)
axis('off')
imshow(im2)
show()
1.2 正向映射和逆向映射
在获得单应矩阵之后，我们的下一个目标是映射图像之间的像素坐标。
操作分为正向映射和逆向映射，归根到底就是上述x’=Hx的一个参数位置的不同。
2. 全景拼接
全景拼接整体流程：
1.根据给定图像集，实现特征匹配
2.通过匹配特征计算图像之间的变换结构
3.利用图像变换结构，实现图像映射
4.针对叠加后的图像，采用APAP之类的算法对齐特征点
5.通过图割算法，自动选取拼接缝
6.根据multi-band blendin策略实现融合
2.1 RANSAC 算法
将上述操作得到的单应性变换图像进行叠加拼接容易出重影现象。这是因为特征点的错误匹配会对通过最小二乘法求解的单应性矩阵带来极大的影响，因此在求解单应性矩阵时，我们需要先通过RANSAC算法实现对于特征点对的筛选，再求解单应性矩阵（将噪声减到最小再求解）。
RANSAC算法流程：
1.随机选择四对匹配特征
2.根据DLT计算单应性矩阵H
3.对所有匹配点，计算映射误差
4.根据误差阈值，确定inliers
5.针对最大inliers集合，重新计算单应性矩阵H
RANSAC算法本质是一个直线拟合算法，我们通过拟合出最佳的直线再设定阈值排除远离拟合直线的噪声点，保留内点。
import numpy as np
import matplotlib.pyplot as plt
import random
import math
# 数据量。
SIZE = 50
# 产生数据。np.linspace 返回一个一维数组，SIZE指定数组长度。
# 数组最小值是0，最大值是10。所有元素间隔相等。
X = np.linspace(0, 10, SIZE)
Y = 3 * X + 10
fig = plt.figure()
# 画图区域分成1行1列。选择第一块区域。
ax1 = fig.add_subplot(1,1, 1)
# 标题
ax1.set_title("RANSAC")
# 让散点图的数据更加随机并且添加一些噪声。
random_x = []
random_y = []
# 添加直线随机噪声
for i in range(SIZE):
random_x.append(X[i] + random.uniform(-0.5, 0.5))
random_y.append(Y[i] + random.uniform(-0.5, 0.5))
# 添加随机噪声
for i in range(SIZE):
random_x.append(random.uniform(0,10))
random_y.append(random.uniform(10,40))
RANDOM_X = np.array(random_x) # 散点图的横轴。
RANDOM_Y = np.array(random_y) # 散点图的纵轴。
# 画散点图。
ax1.scatter(RANDOM_X, RANDOM_Y)
# 横轴名称。
ax1.set_xlabel("x")
# 纵轴名称。
ax1.set_ylabel("y")
# 使用RANSAC算法估算模型
# 迭代最大次数，每次得到更好的估计会优化iters的数值
iters = 100000
# 数据和模型之间可接受的差值
sigma = 0.25
# 最好模型的参数估计和内点数目
best_a = 0
best_b = 0
pretotal = 0
# 希望的得到正确模型的概率
P = 0.99
for i in range(iters):
# 随机在数据中红选出两个点去求解模型
sample_index = random.sample(range(SIZE * 2),2)
x_1 = RANDOM_X[sample_index[0]]
x_2 = RANDOM_X[sample_index[1]]
y_1 = RANDOM_Y[sample_index[0]]
y_2 = RANDOM_Y[sample_index[1]]
# y = ax + b 求解出a，b
a = (y_2 - y_1) / (x_2 - x_1)
b = y_1 - a * x_1
# 算出内点数目
total_inlier = 0
for index in range(SIZE * 2):
y_estimate = a * RANDOM_X[index] + b
if abs(y_estimate - RANDOM_Y[index]) < sigma:
total_inlier = total_inlier + 1
# 判断当前的模型是否比之前估算的模型好
if total_inlier > pretotal:
iters = math.log(1 - P) / math.log(1 - pow(total_inlier / (SIZE * 2), 2))
pretotal = total_inlier
best_a = a
best_b = b
# 判断是否当前模型已经符合超过一半的点
if total_inlier > SIZE:
break
# 用我们得到的最佳估计画图
Y = best_a * RANDOM_X + best_b
# 直线图
ax1.plot(RANDOM_X, Y)
plt.show()
import cv2
import numpy as np
from datetime import datetime
def detectAndCompute(image):
# 计算SIFT特征点和特征向量
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
sift = cv2.xfeatures2d.SIFT_create()
(kps, features) = sift.detectAndCompute(image, None)
kps = np.float32([kp.pt for kp in kps])
return (kps, features)
def matchKeyPoints(kpsA, kpsB, featuresA, featuresB, ratio=0.75, reprojThresh=4.0):
# 匹配
matcher = cv2.BFMatcher()
rawMatches = matcher.knnMatch(featuresA, featuresB, 2)
matches = []
for m in rawMatches:
if len(m) == 2 and m[0].distance < ratio * m[1].distance:
matches.append((m[0].queryIdx, m[0].trainIdx))
# 使用np.float32转化列表
kpsA = np.float32([kpsA[m[0]] for m in matches])
kpsB = np.float32([kpsB[m[1]] for m in matches])
# cv2自带RANSAC筛选状态
(M, status) = cv2.findHomography(kpsA, kpsB, cv2.RANSAC, reprojThresh)
return (M, matches, status)
def drawMatches(imgA, imgB, kpsA, kpsB, matches, status, imageA, imageB):
(hA, wA) = imgA.shape[0:2]
(hB, wB) = imgB.shape[0:2]
# 注意：3通道和uint8类型
drawImg = np.zeros((max(hA, hB), wA + wB, 3), 'uint8')
drawImg[0:hB, 0:wB] = imageB
drawImg[0:hA, wB:] = imageA
for ((queryIdx, trainIdx), s) in zip(matches, status):
if s == 1:
pt1 = (int(kpsB[trainIdx][0]), int(kpsB[trainIdx][1]))
pt2 = (int(kpsA[trainIdx][0]) + wB, int(kpsA[trainIdx][1]))
cv2.line(drawImg, pt1, pt2, (0, 255, 0))
else:
pt1 = (int(kpsB[trainIdx][0]), int(kpsB[trainIdx][1]))
pt2 = (int(kpsA[trainIdx][0]) + wB, int(kpsA[trainIdx][1]))
cv2.line(drawImg, pt1, pt2, (0, 0, 255))
cv2.imwrite("drawMatches.jpg", drawImg)
return
if __name__ == '__main__':
imageA = cv2.imread("test_image/1.jpg")
imageB = cv2.imread("test_image/2.jpg")
(kpsA, featuresA) = detectAndCompute(imageA)
(kpsB, featuresB) = detectAndCompute(imageB)
(M, matches, status) = matchKeyPoints(kpsA, kpsB, featuresA, featuresB)
print(status)
drawMatches(imageA, imageB, kpsA, kpsB, matches, status, imageA, imageB)
2.2 APAP算法
上述通过单应性变换实现的图像映射只能应用于输入图像在一个平面上，或者相机沿着同一个角点旋转拍摄图像，否则容易出现鬼影现象。CVPR2013提出了APAP算法，通过局部单应性变换的思路使两个图像的重叠区域准确对齐。
如下图所示，图上的点表示两张图像之间匹配上的特征点的坐标位置，而红线表示计算出的矩阵。（a）中采用全局的单应性变换操作，会产生很明显的局部偏差；而（c）中采用APAP的方法可以有效的解决局部误差并且推断出全局变换的趋势。
和全局单应性变换相同，求解单应性变换矩阵，但不同的是这里的H是具有位置依赖性的单应性变换矩阵。
这里的N指将图像划分为N个patch，而这里的单应性矩阵需要添加一个权重w（针对每个patch）
权重w的计算由每个像素于该像素所在patch的左上角距离确定，中σ是一个尺度因子，其中距离格子越远的产生的权重越小。这里设计权重的目的我们可以联系上面全局单应性的计算，需要四个特征点的匹配对才能实现投影变换，这时候可能出现patch里没有特征点匹配对的情况，这个时候就需要从其他patch去找特征点匹配对，所以有个权重的概念。而当patch周围的特征点匹配对离它特别远时，权重的概念可以忽略，因此我们采用一个小的值抵消掉权重。
因此局部扭曲变换可以表示为，W是个对角矩阵。
2.3 最小割最大流算法
2.4 multi-band bleing算法
在通过图割算法找到拼接缝后，由于图像光照、噪声等因素，会使拼接缝很突兀。而multi-band bleing算法通过图像融合消除边缘痕迹。
1.对两张图像进行高斯金字塔和拉普拉斯金字塔分解，得到A图和B图的拉普拉斯金字塔。
2.5 全景拼接代码
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
from PCV.localdescriptors import sift
np.seterr(invalid='ignore')
"""
This is the panorama example from section 3.3.
"""
# 设置数据文件夹的路径
featname = ['D:/desktop/cv_task/concate/photo/' + str(i + 1) + '.sift' for i in range(3)]
imname = ['D:/desktop/cv_task/concate/photo/' + str(i + 1) + '.jpg' for i in range(3)]
# 提取特征并匹配使用sift算法
l = {}
d = {}
for i in range(3):
sift.process_image(imname[i], featname[i])
l[i], d[i] = sift.read_features_from_file(featname[i])
matches = {}
for i in range(2):
matches[i] = sift.match(d[i + 1], d[i])
# 可视化匹配
for i in range(2):
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)
# 将匹配转换成齐次坐标点的函数
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
# 估计单应性矩阵
model = homography.RansacModel()
fp, tp = convert_points(1)
H_12 = homography.H_from_ransac(fp, tp, model)[0] # im 1 to 2
fp, tp = convert_points(0)
H_01 = homography.H_from_ransac(fp, tp, model)[0] # im 0 to 1
# 扭曲图像
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)
figure()
imshow(array(im_02, "uint8"))
axis('off')
show()
文章出处登录后可见！
已经登录？立即刷新
