Python计算机视觉-第3章

本章讲解图像之间的变换，以及一些计算变换的实用方法。这些变换可以用于图像 扭曲变形和图像配准。最后，我们将会介绍一个自动创建全景图像的例子。

1、单应性变换

单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。在这里，平 面是指图像或者三维中的平面表面。单应性变换具有很强的实用性，比如图像配准、 图像纠正和纹理扭曲，以及创建全景图像。我们将频繁地使用单应性变换。本质上， 单应性变换 H，按照下面的方程映射二维中的点(齐次坐标意义下):

对于图像平面内(甚至是三维中的点，后面我们会介绍到)的点，齐次坐 标是个非常有用的表示方式。点的齐次坐标是依赖于其尺度定义的，所以， x=[x,y,w]=[αx,αy,αw]=[x/w,y/w,1] 都表示同一个二维点。因此，单应性矩阵 H 也仅 依赖尺度定义，所以，单应性矩阵具有 8 个独立的自由度。我们通常使用 w=1 来归 一化点，这样，点具有唯一的图像坐标 x 和 y。这个额外的坐标使得我们可以简单 地使用一个矩阵来表示变换。

（1）直接线性变换算法

DLT(Direct Linear Transformation，直接线性变换)是给定4个或者更多对应点对 矩阵，来计算单应性矩阵 H 的算法。将单应性矩阵 H 作用在对应点对上，重新写出 该方程，我们可以得到下面的方程:

或者 Ah=0，其中 A 是一个具有对应点对二倍数量行数的矩阵。将这些对应点对方程的系数堆叠到一个矩阵中，我们可以使用 SVD(Singular Value Decomposition， 奇异值分解)算法找到 H 的最小二乘解。下面是该算法的代码。将下面的函数添加 到 homography.py 文件中:

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]

（2）仿射变换

由于仿射变换具有 6 个自由度，因此我们需要三个对应点对来估计矩阵 H。通过将最后两个元素设置为 0，即 h7=h8=0，仿射变换可以用上面的 DLT 算法估计得出。这里我们将使用不同的方法来计算单应性矩阵 H，下面的函数使用对应点对来计算仿射变换矩阵，将其添加到 homograph.py 文件中:

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]

2、图像扭曲

对图像块应用仿射变换，我们将其称为图像扭曲(或者仿射扭曲)。该操作不仅经常应用在计算机图形学中，而且经常出现在计算机视觉算法中。扭曲操作可以使用 SciPy 工具包中的 ndimage 包来简单完成。命令:transformed_im = ndimage.affine_transform(im,A,b,size)

使用如上所示的一个线性变换 A 和一个平移向量 b 来对图像块应用仿射变换。选项 参数 size 可以用来指定输出图像的大小。默认输出图像设置为和原始图像同样大 小。为了研究该函数是如何工作的，我们可以试着运行下面的命令:

from scipy import ndimage
from PIL import Image
from pylab import *

im = array(Image.open('../data/empire.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）图像中的图像

仿射扭曲的一个简单例子是，将图像或者图像的一部分放置在另一幅图像中，使得它们能够和指定的区域或者标记物对齐。代码实现：

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

# example of affine warp of im1 onto im2

im1 = array(Image.open('../data/flickr_image/beatles.png').convert('L'))  #1483 × 957,300x300
im2 = array(Image.open('../data/flickr_image/billboard_for_rent.jpg').convert('L'))  #1475 × 2290,72x72
print(im1.shape[:2])
# set to points  
##tp在该程序中实际上是beatles图片要插在billboard图片上区域的四个顶点坐标(从左上角逆时针依次取点)，tp= （[x1,x2,x3,x4],[y1,y2,y3,y4],[1,1,1,1]])
# tp = array([[264,538,540,264],[40,36,605,605],[1,1,1,1]])
# tp = array([[675,826,826,677],[55,52,281,277],[1,1,1,1]])
tp = array([[580,1400,1400,580],[95,80,1365,1365],[1,1,1,1]])
print(tp)
im3 = warp.image_in_image(im1,im2,tp)
figure()
gray()
subplot(141)
axis('off')
imshow(im1)
subplot(142)
axis('off')
imshow(im2)
subplot(143)
axis('off')
imshow(im3)

tp = array([[1500,1850,1850,1500],[115,110,615,615],[1,1,1,1]])
# set from points to corners of im1
m,n = im1.shape[:2]
fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])
# first triangle
tp2 = tp[:,:3]
fp2 = fp[:,:3]
# compute H
H = homography.Haffine_from_points(tp2,fp2)
im1_t = ndimage.affine_transform(im1,H[:2,:2],
(H[0,2],H[1,2]),im2.shape[:2])
# alpha for triangle
alpha = warp.alpha_for_triangle(tp2,im2.shape[0],im2.shape[1])
im3 = (1-alpha)*im2 + alpha*im1_t
# second triangle
tp2 = tp[:,[0,2,3]]
fp2 = fp[:,[0,2,3]]
# compute H
H = homography.Haffine_from_points(tp2,fp2)
im1_t = ndimage.affine_transform(im1,H[:2,:2],
(H[0,2],H[1,2]),im2.shape[:2])
# alpha for triangle
alpha = warp.alpha_for_triangle(tp2,im2.shape[0],im2.shape[1])
im4 = (1-alpha)*im3 + alpha*im1_t
subplot(144)
imshow(im4)
axis('off')
show()

运行结果：

代码实现：

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

# example of affine warp of im1 onto im2

im1 = array(Image.open('../data/flickr_image/beatles.png').convert('L'))  #1483 × 957,300x300
im2 = array(Image.open('../data/flickr_image/blank_billboard.jpg').convert('L'))  #1475 × 2290,72x72
print(im1.shape[:2])
# set to points  
##tp在该程序中实际上是beatles图片要插在billboard图片上区域的四个顶点坐标(从左上角逆时针依次取点)，tp= （[x1,x2,x3,x4],[y1,y2,y3,y4],[1,1,1,1]])
# tp = array([[264,538,540,264],[40,36,605,605],[1,1,1,1]])
# tp = array([[675,826,826,677],[55,52,281,277],[1,1,1,1]])
tp = array([[580,1415,1165,350],[400,180,2980,2780],[1,1,1,1]])
print(tp)
im3 = warp.image_in_image(im1,im2,tp)
figure()
gray()
subplot(141)
axis('off')
imshow(im1)
subplot(142)
axis('off')
imshow(im2)
subplot(143)
axis('off')
imshow(im3)

# tp = array([[1650,1980,1980,1650],[1515,1510,2015,2015],[1,1,1,1]])
# set from points to corners of im1
m,n = im1.shape[:2]
fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])
# first triangle
tp2 = tp[:,:3]
fp2 = fp[:,:3]
# compute H
H = homography.Haffine_from_points(tp2,fp2)
im1_t = ndimage.affine_transform(im1,H[:2,:2],
(H[0,2],H[1,2]),im2.shape[:2])
# alpha for triangle
alpha = warp.alpha_for_triangle(tp2,im2.shape[0],im2.shape[1])
im3 = (1-alpha)*im2 + alpha*im1_t
# second triangle
tp2 = tp[:,[0,2,3]]
fp2 = fp[:,[0,2,3]]
# compute H
H = homography.Haffine_from_points(tp2,fp2)
im1_t = ndimage.affine_transform(im1,H[:2,:2],
(H[0,2],H[1,2]),im2.shape[:2])
# alpha for triangle
alpha = warp.alpha_for_triangle(tp2,im2.shape[0],im2.shape[1])
im4 = (1-alpha)*im3 + alpha*im1_t
subplot(144)
imshow(im4)
axis('off')
show()

运行结果：

（2）分段仿射扭曲

正如上面的例子所示，三角形图像块的仿射扭曲可以完成角点的精确匹配。让我们 看一下对应点对集合之间最常用的扭曲方式:分段仿射扭曲。给定任意图像的标记 点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，我们可以 将图像和另一幅图像的对应标记点扭曲对应。对于任何图形和图像处理库来说，这 些都是最基本的操作。下面我们来演示一下如何使用 Matplotlib 和 SciPy 来完成该 操作。

一些实例点和三角剖分的代码实现：

 # -*- coding: utf-8 -*-
# import matplotlib.delaunay as md
from scipy.spatial import Delaunay

from pylab import *
from scipy import *

x,y = array(random.standard_normal((2,100)))
# centers,edges,tri,neighbors = md.delaunay(x,y)
tri = Delaunay(np.c_[x,y]).simplices
figure()
gray()
subplot(121)
axis('off')
plot(x,y,'*')
for t in tri:
    t_ext = [t[0], t[1], t[2], t[0]] # add first point to end
    subplot(122)
    plot(x[t_ext],y[t_ext],'r')
plot(x,y,'*')
axis('off')
show()

一些实例点和三角剖分的运行结果：

 

 

使用狄洛克三角剖分标记点进行分段仿射扭曲代码实现：

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

"""
This is the piecewise affine warp example from Section 3.2, Figure 3-5.
"""

# open image to warp
fromim = array(Image.open('../data/sunset_tree.jpg')) #960 × 1280，160 × 160
x, y = meshgrid(range(5), range(6))

x = (fromim.shape[1]/4) * x.flatten()
y = (fromim.shape[0]/5) * y.flatten()

# triangulate
tri = warp.triangulate_points(x, y)

# open image and destination points
im = array(Image.open('../data/turningtorso1.jpg'))#2592 × 3888，240 × 240

figure()
subplot(1, 3, 1)
axis('off')
imshow(im)

tp = loadtxt('../data/turningtorso1_points.txt', 'int')  # destination points
# tp = array([[76,158,158,76],[76,76,158,158],[1,1,1,1]])

# convert points to hom. coordinates (make sure they are of type int)
fp = array(vstack((y, x, ones((1, len(x))))), 'int')
tp = array(vstack((tp[:, 1], tp[:, 0], ones((1, len(tp))))), 'int')

# warp triangles
im = warp.pw_affine(fromim, im, fp, tp, tri)

# plot
subplot(1, 3, 2)
axis('off')
imshow(im)
subplot(1, 3, 3)
axis('off')
imshow(im)
warp.plot_mesh(tp[1], tp[0], tri)

show()

使用狄洛克三角剖分标记点进行分段仿射扭曲运行结果：

（3）图像配准

图像配准是对图像进行变换，使变换后的图像能够在常见的坐标系中对齐。配准可 以是严格配准，也可以是非严格配准。为了能够进行图像对比和更精细的图像分析， 图像配准是一步非常重要的操作。

代码实现：

 # -*- coding: utf-8 -*-
# from PCV.tools import imregistration
import imregistration

"""
This is the face image registration example from Figure 3-6.
Make sure to create a folder 'aligned' under the jkfaces folder.
"""

# load the location of control points
xml_filename = '../data/jkfaces.xml'
points = imregistration.read_points_from_xml(xml_filename)

# register
imregistration.rigid_alignment(points,'../data/jkfaces/')

其中read_points_from_xml函数，rigid_alignment函数和其他函数定义详见imregistration.py：

 # -*- coding: utf-8 -*-
from pylab import *
from xml.dom import minidom
from scipy import linalg,ndimage
from scipy.misc import imsave 
from PIL import Image
import os


def read_points_from_xml(xmlFileName): 
	""" 读取用于人脸对齐的控制点 """

	xmldoc = minidom.parse(xmlFileName)
	facelist = xmldoc.getElementsByTagName('face') 
	faces = {}
	for xmlFace in facelist:
		fileName = xmlFace.attributes['file'].value 
		xf = int(xmlFace.attributes['xf'].value)
		yf = int(xmlFace.attributes['yf'].value)
		xs = int(xmlFace.attributes['xs'].value)
		ys = int(xmlFace.attributes['ys'].value)
		xm = int(xmlFace.attributes['xm'].value)
		ym = int(xmlFace.attributes['ym'].value) 
		faces[fileName] = array([xf, yf, xs, ys, xm, ym])
	return faces


# from scipy import linalg
def compute_rigid_transform(refpoints,points):
	""" 计算用于将点对齐到参考点的旋转、尺度和平移量 """
	A = array([ [points[0], -points[1], 1, 0], 
				[points[1], points[0], 0, 1], 
				[points[2], -points[3], 1, 0], 
				[points[3], points[2], 0, 1],
				[points[4], -points[5], 1, 0], 
				[points[5], points[4], 0, 1]])

	y = array([ refpoints[0], 
				refpoints[1], 
				refpoints[2], 
				refpoints[3], 
				refpoints[4],
				refpoints[5]])

	# 计算最小化 ||Ax-y|| 的最小二乘解
	a,b,tx,ty = linalg.lstsq(A,y)[0]
	R = array([[a, -b], [b, a]]) # 包含尺度的旋转矩阵

	return R,tx,ty

# from PIL import Image
# from scipy import ndimage 
# from scipy.misc import imsave 
# import os
def rigid_alignment(faces,path,plotflag=False): 
	""" 严格对齐图像，并将其保存为新的图像
	path 决定对齐后图像保存的位置 
	设置 plotflag=True，以绘制图像 """

	# 将第一幅图像中的点作为参考点 
	# refpoints = faces.values()[0]  #TypeError: 'dict_values' object is not subscriptable
	refpoints = list(faces.values())[0]
	# 使用仿射变换扭曲每幅图像 
	for face in faces:
		points = faces[face]
	R,tx,ty = compute_rigid_transform(refpoints, points) 
	T = array([[R[1][1], R[1][0]], [R[0][1], R[0][0]]])

	im = array(Image.open(os.path.join(path,face))) 
	im2 = zeros(im.shape, 'uint8')

	# 对每个颜色通道进行扭曲
	for i in range(len(im.shape)):
		im2[:,:,i] = ndimage.affine_transform(im[:,:,i],linalg.inv(T),offset=[-ty,-tx])
	
	if plotflag: 
		imshow(im2) 
		show()

	# 裁剪边界，并保存对齐后的图像 
	h,w = im2.shape[:2]
	border = (w+h)//20

	# 裁剪边界
	# imsave(os.path.join(path, 'aligned/'+face),im2[border:h-border,border:w-border,:])
	imsave(os.path.join(path, 'aligned/'+face),im2[border:h-border,border:w-border,:])

运行结果：

3、创建全景图

在同一位置(即图像的照相机位置相同)拍摄的两幅或者多幅图像是单应性相关的。我们经常使用该约束将很多图像缝补起来，拼成一个大的图像来创建全景图像。在本节中，我们将探讨如何创建全景图像。

（1）RANSAC

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

ransac.py代码实现：

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

## Copyright (c) 2004-2007, Andrew D. Straw. All rights reserved.

## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are
## met:

##     * Redistributions of source code must retain the above copyright
##       notice, this list of conditions and the following disclaimer.

##     * Redistributions in binary form must reproduce the above
##       copyright notice, this list of conditions and the following
##       disclaimer in the documentation and/or other materials provided
##       with the distribution.

##     * Neither the name of the Andrew D. Straw nor the names of its
##       contributors may be used to endorse or promote products derived
##       from this software without specific prior written permission.

## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
## "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
## LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
## A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
## OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
## SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
## LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
## DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
## THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

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
             add point to alsoinliers
    }
    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 = scipy.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:
        # add some outliers
        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 = False
    model = LinearLeastSquaresModel(input_columns,output_columns,debug=debug)

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

    # run RANSAC algorithm
    ransac_fit, ransac_data = ransac(all_data,model,
                                     50, 1000, 7e3, 300, # 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],
                    numpy.dot(A_col0_sorted,linear_fit)[:,0],
                    label='linear fit' )
        pylab.legend()
        pylab.show()

if __name__=='__main__':
    test()
    

ransac.py运行结果：

（2）稳健的单应性矩阵估计

代码实现：

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

import homography, warp
import sift

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

# set paths to data folder
featname = ['../data/Univ'+str(i+1)+'.sift' for i in range(5)] 
imname = ['../data/Univ'+str(i+1)+'.jpg' for i in range(5)]

# extract features and match
l = {}
d = {}
for i in range(5): 
    sift.process_image(imname[i],featname[i])
    l[i],d[i] = sift.read_features_from_file(featname[i])

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

print(matches)

# visualize the matches (Figure 3-11 in the book)
figure()
gray()
for i in range(4):
    subplot(2,2,i+1)
    im1 = array(Image.open(imname[i]))
    im2 = array(Image.open(imname[i+1]))
    sift.plot_matches(im2,im1,l[i+1],l[i],matches[i],show_below=True)

show()

运行结果：

（3）拼接图像

代码实现：

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

import homography, warp
import sift

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

# set paths to data folder
featname = ['../data/Univ'+str(i+1)+'.sift' for i in range(5)] 
imname = ['../data/Univ'+str(i+1)+'.jpg' for i in range(5)]

# extract features and match
l = {}
d = {}
for i in range(5): 
    sift.process_image(imname[i],featname[i])
    l[i],d[i] = sift.read_features_from_file(featname[i])

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

# # visualize the matches (Figure 3-11 in the book)
# figure()
# gray()
# for i in range(4):
#     subplot(2,2,i+1)
#     im1 = array(Image.open(imname[i]))
#     im2 = array(Image.open(imname[i+1]))
#     sift.plot_matches(im2,im1,l[i+1],l[i],matches[i],show_below=True)

# show()


# function to convert the matches to hom. points
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


# estimate the homographies
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    


# warp the images
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()

运行结果：

4、其他

（1）开发环境

         本代码是在mac电脑sublimetxt编辑器python3.7.8下调试出来的，如果你要在windows/linux下编辑器/IDE的其他python版本运行的话，请对代码做相应的调整。

（2）源码和图片

        已经调试过的源码和图片详见：

https://github.com/Abonaventure/pcv-book-code.git   或  https://gitlab.com/Abonaventure/pcv-book-code.git

部分图片可能未上传可在https://pan.baidu.com/share/link?shareid=4059473268&uk=3440544532 或者原书目录http://programmingcomputervision.com下载。