# 图像到图像的映射

## （一）单应性变换

#### 1.1 直接线性变换算法

DLT（Direct Linear Transformation，直接线性变换）是给定４个或者更多对应点对矩阵，来计算单应性矩阵 H 的算法。

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]



#### 1.2 仿射变换

def Haffine_from_points(fp, tp):
""" Find H, affine transformation, such that
tp is affine transf of fp. """

if fp.shape != tp.shape:
raise RuntimeError('number of points do not match')

# condition points
# --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_cond = dot(C1, fp)

# --to points--
m = mean(tp[:2], axis=1)
C2 = C1.copy()  # must use same scaling for both point sets
C2[0][2] = -m[0] / maxstd
C2[1][2] = -m[1] / maxstd
tp_cond = dot(C2, tp)

# conditioned points have mean zero, so translation is zero
A = concatenate((fp_cond[:2], tp_cond[:2]), axis=0)
U, S, V = linalg.svd(A.T)

# create B and C matrices as Hartley-Zisserman (2:nd ed) p 130.
tmp = V[:2].T
B = tmp[:2]
C = tmp[2:4]

tmp2 = concatenate((dot(C, linalg.pinv(B)), zeros((2, 1))), axis=1)
H = vstack((tmp2, [0, 0, 1]))

# decondition
H = dot(linalg.inv(C2), dot(H, C1))

return H / H[2, 2]



## （二）图像扭曲

transformed_im = ndimage.affine_transform(im,A,b,size)

# -*- coding: utf-8 -*-
from scipy import ndimage
from PIL import Image
from pylab import *

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)

im = array(Image.open('D:\\Python\\chapter3\\jimei.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)
title(u'(a)原始图像', fontproperties=font)
axis('off')
imshow(im)
subplot(122)
title(u'(b)使用 ndimage.affine_transform()函数扭曲后的图像', fontproperties=font)
axis('off')
imshow(im2)
show()


#### 2.1 图像中的图像

# -*- coding: utf-8 -*-
from PCV.geometry import warp, homography
from PIL import Image
from pylab import *
from scipy import ndimage

# 将im1仿射扭曲到im2的指定位置
im1 = array(Image.open('D:\\Python\\chapter3\\star.jpg').convert('L'))
# 选定指定的位置，即目标点
tp = array([[264,538,540,264],[40,36,605,605],[1,1,1,1]])
# 调用的warp.py的image_in_image函数，从而实现仿射变换
im3 = warp.image_in_image(im1, im2, tp)
figure()
gray()
subplot(131)
axis('off')
imshow(im1)
subplot(132)
axis('off')
imshow(im2)
subplot(133)
axis('off')
imshow(im3)
axis('off')
show()


#### 2.2 图像配准

f为二维空间坐标变换（如仿射变换），g为一维亮度或其他度量值变换。

1. 特征提取
2. 特征匹配
3. 估计变换模型
4. 图像重采样及变换

（1）特征提取：

a. Harris算法

b. Susan算法
susan算法使用一个圆形的模板在图像上滑动，将位于圆形模板中心的待检测的像素点称为核心点。假设图像为非纹理，核心点领域被划分为两个区域：其一为亮度值等于（或相似于）核心点亮度的区域，称为核值相似区，其二为亮度值不相似于核心点亮度的区域。

c. Harris-Laplace算法
Harris算子能最稳定地在图像旋转、光照变化、透视变换条件下提取二维平面特征点，但在三维尺度空间中，Harris探测子的重复探测性能不好，不同尺度Harris特征点存在位置误差，Harris探测子不具有尺度和仿射不变性。而三维尺度空间中最稳定高效的特征尺度探测算子是归一化的Laplace算子。Harris-Laplace的特征点具有尺度和旋转不变得特性，且对光照变换和小范围视角变化具有稳定性。

d. SIFT特征点提取

e. SURF特征点提取

（ 2 ) 特征匹配

a. 对特征作描述

b. 利用相似度准则进行特征匹配

SIFT特征描述子:

SURF特征描述子:

DAISY特征描述子:

（3）估计变换模型

a.刚体变换模型

b.仿射变换模型

c.投影变换模型

d.非线性变换模型

F表示把第一幅图像映射到第二幅图像上的任意一种函数形式。典型的非线性变换如多项式变换，在2D空间中，多项式函数可写成如下形式：

（4）图像重采样及变换

（1） 基于点：

a. SIFT算法
SIFT特征匹配算法包括两个阶段:SIFT特征的生成与SIFT特征向量的匹配。
SIFT特征向量的生成算法包括四步：
1.尺度空间极值检测，以初步确定关键点位置和所在尺度。
2.拟和三维二次函数精确确定位置和尺度，同时去除低对比度的关键点和不稳定的边缘响应点。
3.利用关键点领域像素的梯度方向分布特性为每个关键点指定参数方向，使算子具备旋转不变性。
4.生成SIFT特征向量。
SIFT特征向量的匹配

b.ASIFT算法

c.SUFR算法

1.计算原图像的积分图像
2.用不同尺寸的框状滤波器来计算不同阶以及不同层上的每个点图像点的行列式。一般计算4阶4层
3.在3维(x,y,S）尺度空间中，在每个3x3x3的局部区域里，进行非最大值抑制。只有比邻近的26个点的响应值都大的点才被选为兴趣点。

（2）基于边缘：

1.根据边缘的相似度
2.提取边缘上的控制点，如曲率比较大的点等，然后用这些点来进行匹配。
3.将边缘拟合成直线，然后匹配直线。

# -*- coding: utf-8 -*-
import numpy as np
import cv2

def sift_kp(image):
gray_image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
sift = cv2.xfeatures2d_SIFT.create()
kp, des = sift.detectAndCompute(image, None)
kp_image = cv2.drawKeypoints(gray_image, kp, None)
return kp_image, kp, des

def get_good_match(des1, des2):
bf = cv2.BFMatcher()
matches = bf.knnMatch(des1, des2, k=2)
good = []
for m, n in matches:
if m.distance < 0.75 * n.distance:
good.append(m)
return good

def siftImageAlignment(img1, img2):
_, kp1, des1 = sift_kp(img1)
_, kp2, des2 = sift_kp(img2)
goodMatch = get_good_match(des1, des2)
if len(goodMatch) > 4:
ptsA = np.float32([kp1[m.queryIdx].pt for m in goodMatch]).reshape(-1, 1, 2)
ptsB = np.float32([kp2[m.trainIdx].pt for m in goodMatch]).reshape(-1, 1, 2)
ransacReprojThreshold = 4
H, status = cv2.findHomography(ptsA, ptsB, cv2.RANSAC, ransacReprojThreshold);
# 其中H为求得的单应性矩阵矩阵
# status则返回一个列表来表征匹配成功的特征点。
# ptsA,ptsB为关键点
# cv2.RANSAC, ransacReprojThreshold这两个参数与RANSAC有关
imgOut = cv2.warpPerspective(img2, H, (img1.shape[1], img1.shape[0]),
flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP)
return imgOut, H, status

while img1.shape[0] > 1000 or img1.shape[1] > 1000:
img1 = cv2.resize(img1, None, fx=0.5, fy=0.5, interpolation=cv2.INTER_AREA)
while img2.shape[0] > 1000 or img2.shape[1] > 1000:
img2 = cv2.resize(img2, None, fx=0.5, fy=0.5, interpolation=cv2.INTER_AREA)

result, _, _ = siftImageAlignment(img1, img2)
allImg = np.concatenate((img1, img2, result), axis=1)
cv2.namedWindow('1', cv2.WINDOW_NORMAL)
cv2.namedWindow('2', cv2.WINDOW_NORMAL)
cv2.namedWindow('Result', cv2.WINDOW_NORMAL)
cv2.imshow('1', img1)
cv2.imshow('2', img2)
cv2.imshow('Result', result)
# cv2.imshow('Result',allImg)
if cv2.waitKey(200000) & 0xff == ord('q'):
cv2.destroyAllWindows()
cv2.waitKey(1)


# Python 2/3 compatibility
from __future__ import print_function

import numpy as np
import cv2 as cv

# built-in modules
import itertools as it

# local modules
from common import Timer
from find_obj import init_feature, filter_matches, explore_match

'''

Ai - is an affine transform matrix from skew_img to img
'''
h, w = img.shape[:2]
A = np.float32([[1, 0, 0], [0, 1, 0]])
if phi != 0.0:
s, c = np.sin(phi), np.cos(phi)
A = np.float32([[c,-s], [ s, c]])
corners = [[0, 0], [w, 0], [w, h], [0, h]]
tcorners = np.int32( np.dot(corners, A.T) )
x, y, w, h = cv.boundingRect(tcorners.reshape(1,-1,2))
A = np.hstack([A, [[-x], [-y]]])
img = cv.warpAffine(img, A, (w, h), flags=cv.INTER_LINEAR, borderMode=cv.BORDER_REPLICATE)
if tilt != 1.0:
s = 0.8*np.sqrt(tilt*tilt-1)
img = cv.GaussianBlur(img, (0, 0), sigmaX=s, sigmaY=0.01)
img = cv.resize(img, (0, 0), fx=1.0/tilt, fy=1.0, interpolation=cv.INTER_NEAREST)
A[0] /= tilt
if phi != 0.0 or tilt != 1.0:
h, w = img.shape[:2]
Ai = cv.invertAffineTransform(A)

'''
affine_detect(detector, img, mask=None, pool=None) -> keypoints, descrs

Apply a set of affine transformations to the image, detect keypoints and
reproject them into initial image coordinates.
See http://www.ipol.im/pub/algo/my_affine_sift/ for the details.

ThreadPool object may be passed to speedup the computation.
'''
params = [(1.0, 0.0)]
for t in 2**(0.5*np.arange(1,6)):
for phi in np.arange(0, 180, 72.0 / t):
params.append((t, phi))

def f(p):
t, phi = p
timg, tmask, Ai = affine_skew(t, phi, img)
for kp in keypoints:
x, y = kp.pt
kp.pt = tuple( np.dot(Ai, (x, y, 1)) )
if descrs is None:
descrs = []
return keypoints, descrs

keypoints, descrs = [], []
if pool is None:
ires = it.imap(f, params)
else:
ires = pool.imap(f, params)

for i, (k, d) in enumerate(ires):
print('affine sampling: %d / %d\r' % (i+1, len(params)), end='')
keypoints.extend(k)
descrs.extend(d)

print()
return keypoints, np.array(descrs)

if __name__ == '__main__':
print(__doc__)

import sys, getopt
opts, args = getopt.getopt(sys.argv[1:], '', ['feature='])
opts = dict(opts)
'''
--feature:  sift/surf/orb/akaze
'''
feature_name = opts.get('--feature', 'brisk-flann')
try:
fn1, fn2 = args
except:
fn1 = 'D:\\Python\\chapter3\\building.jpg'
fn2 = 'D:\\Python\\chapter3\\building2.jpg'

detector, matcher = init_feature(feature_name)

if img1 is None:
sys.exit(1)

if img2 is None:
sys.exit(1)

if detector is None:
print('unknown feature:', feature_name)
sys.exit(1)
print('using', feature_name)

kp1, desc1 = affine_detect(detector, img1, pool=pool)
kp2, desc2 = affine_detect(detector, img2, pool=pool)
print('img1 - %d features, img2 - %d features' % (len(kp1), len(kp2)))

def match_and_draw(win):
with Timer('matching'):
raw_matches = matcher.knnMatch(desc1, trainDescriptors=desc2, k=2) #2
p1, p2, kp_pairs = filter_matches(kp1, kp2, raw_matches)
if len(p1) >= 4:
H, status = cv.findHomography(p1, p2, cv.RANSAC, 5.0)
print('%d / %d  inliers/matched' % (np.sum(status), len(status)))
# do not draw outliers (there will be a lot of them)
kp_pairs = [kpp for kpp, flag in zip(kp_pairs, status) if flag]
else:
H, status = None, None
print('%d matches found, not enough for homography estimation' % len(p1))
explore_match(win, img1, img2, kp_pairs, None, H)

match_and_draw('affine find_obj')
cv.waitKey()
cv.destroyAllWindows()



## （三）创建全景图

#### 3.1 RANSAC

RANSAC 是“RANdom SAmple Consensus”（随机一致性采样）的缩写。该方法是用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型，例如点集之间的单应性矩阵，RANSAC 基本的思想是，数据中包含正确的点和噪声点，合理的模型应该能够在描述正确数据点的同时摒弃噪声点。

RANSAC的基本假设是：
（1）数据由“局内点”组成，例如：数据的分布可以用一些模型参数来解释；
（2）“局外点”是不能适应该模型的数据；
（3）除此之外的数据属于噪声。

RANSAC也做了以下假设：给定一组（通常很小的）局内点，存在一个可以估计模型参数的过程；而该模型能够解释或者适用于局内点。

RANSAC算法的输入是一组观测数据，一个可以解释或者适应于观测数据的参数化模型，一些可信的参数。
RANSAC通过反复选择数据中的一组随机子集来达成目标。被选取的子集被假设为局内点，并用下述方法进行验证：

1. 首先我们先随机假设一小组局内点为初始值。然后用此局内点拟合一个模型，此模型适应于假设的局内点，所有的未知参数都能从假设的局内点计算得出。
2. 用1中得到的模型去测试所有的其它数据，如果某个点适用于估计的模型，认为它也是局内点，将局内点扩充。
3. 如果有足够多的点被归类为假设的局内点，那么估计的模型就足够合理。
4. 然后，用所有假设的局内点去重新估计模型，因为此模型仅仅是在初始的假设的局内点估计的，后续有扩充后，需要更新。
5. 最后，通过估计局内点与模型的错误率来评估模型。

1、要么因为局内点太少，还不如上一次的模型，而被舍弃，
2、要么因为比现有的模型更好而被选用。

RANSAC原理
RANSAC目的是找到最优的参数矩阵使得满足该矩阵的数据点个数最多，通常令h33=1来归一化矩阵。由于单应性矩阵有8个未知参数，至少需要8个线性方程求解，对应到点位置信息上，一组点对可以列出两个方程，则至少包含4组匹配点对。

RANSAC算法从匹配数据集中随机抽出4个样本并保证这4个样本之间不共线，计算出单应性矩阵，然后利用这个模型测试所有数据，并计算满足这个模型数据点的个数与投影误差(即代价函数)，若此模型为最优模型，则对应的代价函数最小。

RANSAC 的标准例子：用一条直线拟合带有噪声数据的点集。简单的最小二乘在该 例子中可能会失效，但是 RANSAC 能够挑选出正确的点，然后获取能够正确拟合的直线。

import numpy
import scipy  # use numpy if scipy unavailable
import scipy.linalg  # use numpy if scipy unavailable

def ransac(data, model, n, k, t, d, debug=False, return_all=False):
"""fit model parameters to data using the RANSAC algorithm

This implementation written from pseudocode found at
http://en.wikipedia.org/w/index.php?title=RANSAC&oldid=116358182

{{{
Given:
data - a set of observed data points
model - a model that can be fitted to data points
n - the minimum number of data values required to fit the model
k - the maximum number of iterations allowed in the algorithm
t - a threshold value for determining when a data point fits a model
d - the number of close data values required to assert that a model fits well to data
Return:
bestfit - model parameters which best fit the data (or nil if no good model is found)
iterations = 0
bestfit = nil
besterr = something really large
while iterations < k {
maybeinliers = n randomly selected values from data
maybemodel = model parameters fitted to maybeinliers
alsoinliers = empty set
for every point in data not in maybeinliers {
if point fits maybemodel with an error smaller than t
}
if the number of elements in alsoinliers is > d {
% this implies that we may have found a good model
% now test how good it is
bettermodel = model parameters fitted to all points in maybeinliers and alsoinliers
thiserr = a measure of how well model fits these points
if thiserr < besterr {
bestfit = bettermodel
besterr = thiserr
}
}
increment iterations
}
return bestfit
}}}
"""
iterations = 0
bestfit = None
besterr = numpy.inf
best_inlier_idxs = None
while iterations < k:
maybe_idxs, test_idxs = random_partition(n, data.shape[0])
maybeinliers = data[maybe_idxs, :]
test_points = data[test_idxs]
maybemodel = model.fit(maybeinliers)
test_err = model.get_error(test_points, maybemodel)
also_idxs = test_idxs[test_err < t]  # select indices of rows with accepted points
alsoinliers = data[also_idxs, :]
if debug:
print 'test_err.min()', test_err.min()
print 'test_err.max()', test_err.max()
print 'numpy.mean(test_err)', numpy.mean(test_err)
print 'iteration %d:len(alsoinliers) = %d' % (
iterations, len(alsoinliers))
if len(alsoinliers) > d:
betterdata = numpy.concatenate((maybeinliers, alsoinliers))
bettermodel = model.fit(betterdata)
better_errs = model.get_error(betterdata, bettermodel)
thiserr = numpy.mean(better_errs)
if thiserr < besterr:
bestfit = bettermodel
besterr = thiserr
best_inlier_idxs = numpy.concatenate((maybe_idxs, also_idxs))
iterations += 1
if bestfit is None:
raise ValueError("did not meet fit acceptance criteria")
if return_all:
return bestfit, {'inliers': best_inlier_idxs}
else:
return bestfit

def random_partition(n, n_data):
"""return n random rows of data (and also the other len(data)-n rows)"""
all_idxs = numpy.arange(n_data)
numpy.random.shuffle(all_idxs)
idxs1 = all_idxs[:n]
idxs2 = all_idxs[n:]
return idxs1, idxs2

class LinearLeastSquaresModel:
"""linear system solved using linear least squares

This class serves as an example that fulfills the model interface
needed by the ransac() function.

"""

def __init__(self, input_columns, output_columns, debug=False):
self.input_columns = input_columns
self.output_columns = output_columns
self.debug = debug

def fit(self, data):
A = numpy.vstack([data[:, i] for i in self.input_columns]).T
B = numpy.vstack([data[:, i] for i in self.output_columns]).T
x, resids, rank, s = numpy.linalg.lstsq(A, B)
return x

def get_error(self, data, model):
A = numpy.vstack([data[:, i] for i in self.input_columns]).T
B = numpy.vstack([data[:, i] for i in self.output_columns]).T
B_fit = scipy.dot(A, model)
err_per_point = numpy.sum((B - B_fit) ** 2, axis=1)  # sum squared error per row
return err_per_point

def test():
# generate perfect input data

n_samples = 500
n_inputs = 1
n_outputs = 1
A_exact = 20 * numpy.random.random((n_samples, n_inputs))
perfect_fit = 60 * numpy.random.normal(size=(n_inputs, n_outputs))  # the model
B_exact = scipy.dot(A_exact, perfect_fit)
assert B_exact.shape == (n_samples, n_outputs)

# add a little gaussian noise (linear least squares alone should handle this well)
A_noisy = A_exact + numpy.random.normal(size=A_exact.shape)
B_noisy = B_exact + numpy.random.normal(size=B_exact.shape)

if 1:
n_outliers = 100
all_idxs = numpy.arange(A_noisy.shape[0])
numpy.random.shuffle(all_idxs)
outlier_idxs = all_idxs[:n_outliers]
non_outlier_idxs = all_idxs[n_outliers:]
A_noisy[outlier_idxs] = 20 * numpy.random.random((n_outliers, n_inputs))
B_noisy[outlier_idxs] = 50 * numpy.random.normal(size=(n_outliers, n_outputs))

# setup model

all_data = numpy.hstack((A_noisy, B_noisy))
input_columns = range(n_inputs)  # the first columns of the array
output_columns = [n_inputs + i for i in range(n_outputs)]  # the last columns of the array
debug = True
model = LinearLeastSquaresModel(input_columns, output_columns, debug=debug)

linear_fit, resids, rank, s = numpy.linalg.lstsq(all_data[:, input_columns], all_data[:, output_columns])

# run RANSAC algorithm
ransac_fit, ransac_data = ransac(all_data, model,
5, 5000, 7e4, 50,  # misc. parameters
debug=debug, return_all=True)
if 1:
import pylab

sort_idxs = numpy.argsort(A_exact[:, 0])
A_col0_sorted = A_exact[sort_idxs]  # maintain as rank-2 array

if 1:
pylab.plot(A_noisy[:, 0], B_noisy[:, 0], 'k.', label='data')
pylab.plot(A_noisy[ransac_data['inliers'], 0], B_noisy[ransac_data['inliers'], 0], 'bx',
label='RANSAC data')
else:
pylab.plot(A_noisy[non_outlier_idxs, 0], B_noisy[non_outlier_idxs, 0], 'k.', label='noisy data')
pylab.plot(A_noisy[outlier_idxs, 0], B_noisy[outlier_idxs, 0], 'r.', label='outlier data')
pylab.plot(A_col0_sorted[:, 0],
numpy.dot(A_col0_sorted, ransac_fit)[:, 0],
label='RANSAC fit')
pylab.plot(A_col0_sorted[:, 0],
numpy.dot(A_col0_sorted, perfect_fit)[:, 0],
label='exact system')
pylab.plot(A_col0_sorted[:, 0],
nu  mpy.dot(A_col0_sorted, linear_fit)[:, 0],
label='linear fit')
pylab.legend()
pylab.show()
if __name__ == '__main__':
test()


RANSAC算法的输入为

1. 观测数据 （包括内群与外群的数据）
2. 符合部分观测数据的模型 （与内群相符的模型）
3. 最少符合模型的内群数量
4. 判断数据是否符合模型的阈值 （数据与模型之间的误差容忍度）
5. 迭代运算次数 （抽取多少次随机内群）

1. 最符合数据的模型参数 （如果内群数量小于输入第三条则判断为数据不存在此模型）
2. 内群集 （符合模型的数据）

#### 3.2 稳健的单应性矩阵估计

featname = ['D:\\Python\\chapter3\\building\\building' + str(i + 1) + '.sift' for i in range(5)]
imname = ['D:\\Python\\chapter3\\building\\building' + str(i + 1) + '.jpg' for i in range(5)]

# 提取特征并匹配使用sift算法
l = {}
d = {}
for i in range(5):
sift.process_image(imname[i], featname[i])

matches = {}
for i in range(4):
matches[i] = sift.match(d[i + 1], d[i])


#### 3.3 拼接图像

1. 特征点匹配：找到素材图片中共有的图像部分。
2. 图片匹配：连接匹配的特征点，估算图像间几何方面的变换。

1. 全景图矫直：矫正拍摄图片时相机的相对3D旋转，主要原因是拍摄图片时相机很可能并不在同一水平线上，并且存在不同程度的倾斜，略过这一步可能导致全景图变成波浪形状。
2. 图像均衡补偿：全局平衡所有图片的光照和色调。
3. 图像频段融合：步骤4之后仍然会存在图像之间衔接边缘、晕影效果（图像的外围部分的亮度或饱和度比中心区域低）、视差效果（因为相机透镜移动导致）。

def panorama(H, fromim, toim, padding=2400, delta=2400):
""" Create horizontal panorama by blending two images
using a homography H (preferably estimated using RANSAC).
The result is an image with the same height as toim. 'padding'
specifies number of fill pixels and 'delta' additional translation. """

# check if images are grayscale or color
is_color = len(fromim.shape) == 3

# homography transformation for geometric_transform()
def transf(p):
p2 = dot(H, [p[0], p[1], 1])
return (p2[0] / p2[2], p2[1] / p2[2])

if H[1, 2] < 0:  # fromim is to the right
print ('warp - right')
# transform fromim
if is_color:
# pad the destination image with zeros to the right
toim_t = hstack((toim, zeros((toim.shape[0], padding, 3))))
fromim_t = zeros((toim.shape[0], toim.shape[1] + padding, toim.shape[2]))
for col in range(3):
fromim_t[:, :, col] = ndimage.geometric_transform(fromim[:, :, col],
else:
# pad the destination image with zeros to the right
fromim_t = ndimage.geometric_transform(fromim, transf,
else:
print ('warp - left')
H_delta = array([[1, 0, 0], [0, 1, -delta], [0, 0, 1]])
H = dot(H, H_delta)
# transform fromim
if is_color:
# pad the destination image with zeros to the left
toim_t = hstack((zeros((toim.shape[0], padding, 3)), toim))
fromim_t = zeros((toim.shape[0], toim.shape[1] + padding, toim.shape[2]))
for col in range(3):
fromim_t[:, :, col] = ndimage.geometric_transform(fromim[:, :, col],
else:
# pad the destination image with zeros to the left
fromim_t = ndimage.geometric_transform(fromim,

# blend and return (put fromim above toim)
if is_color:
# all non black pixels
alpha = ((fromim_t[:, :, 0] * fromim_t[:, :, 1] * fromim_t[:, :, 2]) > 0)
for col in range(3):
toim_t[:, :, col] = fromim_t[:, :, col] * alpha + toim_t[:, :, col] * (1 - alpha)
else:
alpha = (fromim_t > 0)
toim_t = fromim_t * alpha + toim_t * (1 - alpha)


# -*- coding: utf-8 -*-
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

"""
This is the panorama example from section 3.3.
"""

# 设置数据文件夹的路径
featname = ['D:\\Python\\chapter3\\building\\building' + str(i + 1) + '.sift' for i in range(5)]
imname = ['D:\\Python\\chapter3\\building\\building' + str(i + 1) + '.jpg' for i in range(5)]

# 提取特征并匹配使用sift算法
l = {}
d = {}
for i in range(5):
sift.process_image(imname[i], featname[i])

matches = {}
for i in range(4):
matches[i] = sift.match(d[i + 1], d[i])

# 可视化匹配
for i in range(4):
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

tp, fp = convert_points(2)  # NB: reverse order
H_32 = homography.H_from_ransac(fp, tp, model)[0]  # im 3 to 2

tp, fp = convert_points(3)  # NB: reverse order
H_43 = homography.H_from_ransac(fp, tp, model)[0]  # im 4 to 3

# 扭曲图像
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)

im1 = array(Image.open(imname[3]), "f")
im_32 = warp.panorama(H_32,im1,im_02,delta,delta)

im1 = array(Image.open(imname[4]), "f")
im_42 = warp.panorama(dot(H_32,H_43),im1,im_32,delta,2*delta)

figure()
imshow(array(im_42, "uint8"))
axis('off')
show()



©️2019 CSDN 皮肤主题: 技术工厂 设计师: CSDN官方博客