用随机采样一致性方法求单应性矩阵H

先看结果：

需要配准的图片为img1 和 img2，如下图：

需要配准，首先需要找关键点：

注意：这里显示的关键点，是配准后，发现可以找到匹配点的关键点，并没有显示原始关键点。

然后用随机一致性(RANRAC)方法，找到这些关键点最合适的配对，如下图：

配准后，得到单应性矩阵H：

然后使用单应性矩阵对img1进行变换：

变换后，发现原img1上的海面部分不见了，为了将海面部分也显示出来，所以，我们先将img1右移640个像素，然后再通过H变换，得：

将变换后的图像与img2相加，就会发现这两张图片已经配准了。

因为img1是先右移640像素，才进行H变换的。

所以，这里，我们也要将img2右移640像素。

下面我们看看详细的原理：

Affine invariant feature-based image matching sample.
This sample is similar to find_obj.py, but uses the affine transformation
space sampling technique, called ASIFT [1]. While the original implementation
is based on SIFT, you can try to use SURF or ORB detectors instead. Homography RANSAC
is used to reject outliers. Threading is used for faster affine sampling.

[1] http://www.ipol.im/pub/algo/my_affine_sift/

USAGE
  asift.py [--feature=<sift|surf|orb|brisk>[-flann]] [ <image1> <image2> ]

  --feature  - Feature to use. Can be sift, surf, orb or brisk. Append '-flann'
               to feature name to use Flann-based matcher instead bruteforce.

  Press left mouse button on a feature point to see its matching point.

affine transformation space sampling technique:刚体变换空间采样技术

Homography RANSAC 用来拒绝例外点。

所以，这里面涉及到两个技术：一个是SIFT,一个是Homography RANSAC

要说清楚Homography RANSAC，就需要先了解RANSAC

这个思路其实是估计出了k个模型，在这k个模型中，有95%的概率，存在一个模型符合所有的正确数据（局内点）。

这里解释一下这个局内点与局外点。

局外点就是噪声数据。

局内点就是正常数据。

所以，RANSAC方法通过概率的形式，去除了噪声点对模型的影响。并且能够判断出来，哪些是局外点。

说清楚了RANSAC，我们再看Homography RANSAC

Homography：单应性。这是个啥？

这里：单应，单独对应，一一对应，点与点是一一对应的关系。

图中的H就是单应行矩阵。而用单应性矩阵乘以世界坐标系下点的坐标（X,Y,Z)就会得到像素坐标系下的坐标(u,v),而这就是单应行变换。

我们看看单应性变换能做的事情：

我们做图像匹配的时候，就是找两幅图上像素点之间的单应行矩阵。

找到这个，两幅图的匹配就搞定了。

匹配的过程：

1 找到a,b图上的关键点集合A，B。

2 求集合A和集合B的单应性矩阵。

3 利用单应行矩阵，将b图转化，然后叠加在A图上。

单应性矩阵求解

这显然有9个未知量。

因为我们有集合A，B是已知的，

我们可以利用这些点，建立方程组，然后求解方程组的方式求取得到这9个值。

也可以用矩阵的计算求取得到单应性矩阵的值。

也可以通过随机采样一致性(RANSC)方法进行竞选，选出一个最好的模型。

参考：
https://blog.csdn.net/qq_45467083/article/details/105457198?utm_medium=distribute.pc_relevant.none-task-blog-2~default~baidujs_title~default-0.control&spm=1001.2101.3001.4242

参考：
https://blog.csdn.net/lhanchao/article/details/52849446

具体实现代码如下：

#!/usr/bin/env python

'''
Affine invariant feature-based image matching sample.

This sample is similar to find_obj.py, but uses the affine transformation
space sampling technique, called ASIFT [1]. While the original implementation
is based on SIFT, you can try to use SURF or ORB detectors instead. Homography RANSAC
is used to reject outliers. Threading is used for faster affine sampling.

[1] http://www.ipol.im/pub/algo/my_affine_sift/

USAGE
asift.py [--feature=<sift|surf|orb|brisk>[-flann]] [ <image1> <image2> ]

--feature - Feature to use. Can be sift, surf, orb or brisk. Append '-flann'
to feature name to use Flann-based matcher instead bruteforce.

Press left mouse button on a feature point to see its matching point.
'''

# Python 2/3 compatibility
from __future__ import print_function

import numpy as np
import cv2 as cv

# built-in modules
import itertools as it
from multiprocessing.pool import ThreadPool

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


def affine_skew(tilt, phi, img, mask=None):
'''
affine_skew(tilt, phi, img, mask=None) -> skew_img, skew_mask, Ai

Ai - is an affine transform matrix from skew_img to img
'''
h, w = img.shape[:2]
if mask is None:
mask = np.zeros((h, w), np.uint8)
mask[:] = 255
A = np.float32([[1, 0, 0], [0, 1, 0]])
if phi != 0.0:
phi = np.deg2rad(phi)
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]
mask = cv.warpAffine(mask, A, (w, h), flags=cv.INTER_NEAREST)
Ai = cv.invertAffineTransform(A)
return img, mask, Ai


def affine_detect(detector, img, mask=None, pool=None):
'''
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)
keypoints, descrs = detector.detectAndCompute(timg, tmask)
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)


def main():
import sys, getopt
opts, args = getopt.getopt(sys.argv[1:], '', ['feature='])
opts = dict(opts)
feature_name = opts.get('--feature', 'brisk-flann')
try:
fn1, fn2 = args
except:
fn1 = 'aero3.jpg'
fn2 = 'aero1.jpg'

img1 = cv.imread(cv.samples.findFile(fn1), cv.IMREAD_GRAYSCALE)
img2 = cv.imread(cv.samples.findFile(fn2), cv.IMREAD_GRAYSCALE)
detector, matcher = init_feature(feature_name)

if img1 is None:
print('Failed to load fn1:', fn1)
sys.exit(1)

if img2 is None:
print('Failed to load fn2:', fn2)
sys.exit(1)

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

print('using', feature_name)

pool=ThreadPool(processes = cv.getNumberOfCPUs())
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)
#Hi = cv.invert(H)
print(H)
# 为了适应图片1变换后的图像大小，我们将原图平移，并放到更大的画布上
H2 = np.array([[1.,0.,640],[0.,1.,0.],[0.,0.,1.]])
img1H=cv.warpPerspective(img1,H,[img1.shape[1],img1.shape[0]])

# 将img1平移并利用单应性矩阵进行变换
img1H2=cv.warpPerspective(img1,np.dot(H2,H),[640+2*img1.shape[1],img1.shape[0]])
# 将图片2平移
img2H2=cv.warpPerspective(img2,H2,[640+2*img2.shape[1],img2.shape[0]])
#img2H=cv.warpPerspective(img2,Hi,[img1.shape[1],img1.shape[0]])
cv.imshow('img1H',np.vstack([img1,img1H,img2]))
cv.waitKey()
cv.imshow('img1H2',np.vstack([img1H2]))
cv.waitKey()
cv.imshow('img2H2',np.vstack([img2H2]))
cv.waitKey()
cv.imshow('img-addweight',cv.addWeighted(img1H2,0.5,img2H2,0.5,0))
cv.waitKey()
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()
print('Done')


if __name__ == '__main__':
print(__doc__)
main()
cv.destroyAllWindows()