计算机视觉编程 第三章

前言

本专栏按《python计算机视觉编程 ——Jan Erik Solem》一书为参考，第三章介绍图像间的变换，如图像扭曲、图像位置映射、图像匹配和构建全景图等。本质上是以线性代数变换为基础的图像映射方法

3.1 单应性变换

单应性变换H描述了在平面上的两个二维坐标系之间的映射关系，能够用于实现图像的投影、变形、配准以及多视角图像合成等操作。它可以被表示为一个3x3的矩阵，以如下形式进行平面间映射：$$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right]=\left[\begin{array}{c}{h}_1{h}_2{h}_3\\ h_4{h}_5{h}_6\\ h_7{h}_8{h}_9\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right]或{\mathbf{x}}^{\prime }=\mathbf{Hx}$$由于点的齐次坐标依赖于其尺度定义，因此同一个二维点$\mathbf{x}$有$[x,y,w]=[\alpha x,\alpha y,\alpha w]=[x/w,y/w,1]$即八个独立自由度，单应性变换H也有相同性质，通常使用$w=1$来归一化，这样$n$个二维点就转化为齐次坐标下的$3\times n$数组。当然不同的特征处理也会要求不同的数据形式，下面是进行归一化和转换齐次坐标的代码

def normalize(points):
    """ Normalize a collection of points in 
        homogeneous coordinates so that last row = 1. """

    for row in points:
        row /= points[-1]
    return points
    
    
def make_homog(points):
    """ Convert a set of points (dim*n array) to 
        homogeneous coordinates. """
        
    return vstack((points,ones((1,points.shape[1])))) 

对于本章主要介绍三种投影变换：

• 仿射变换：$$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{c}{a}_1{a}_2{t}_x\\ a_3{a}_4{t}_y\\ 001\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]或{\mathbf{x}}^{\prime }=\left[\begin{array}{c}At\\ 01\end{array}\right]x$$它有可逆矩阵$A$和平移向量$t=[t_x,t_y]$组成，可以包括平移、旋转、缩放和剪切等操作，有6个自由度
• 相似变换：$$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{ccc}s\cos (\theta )& -s\sin (\theta )& t_x\\ s\sin (\theta )& s\cos (\theta )& t_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]或{\mathbf{x}}^{\prime }=\left[\begin{array}{c}sRt\\ 01\end{array}\right]x$$是包含尺度变化的二维刚题变换，向量$s$指定变换的尺度，若它为1则为刚体变换。$R$是角度为$\theta$的旋转矩阵。相似变换是仿射变换的特例，保持了点的比例关系和角度。在相似变换中，平移、旋转和缩放仍然是主要操作，但没有剪切
• 完全投影变换：公式同单应性变换，是一种更通用的变换，可以包含透视效果

3.1.1 直接线性变换算法

一个完全射影变换有8个自由度，根据对应点约束，每个对应点对应$x$和$y$两个方程，因此计算单应性矩阵$H$需要四个对应点对。直接线性变换算法DLT就是给定四个或者多个对应点对矩阵计算$H$ 的算法。将单应性矩阵$H$ 作用在对应点对上就得到了一个具有对应点对二倍数量行数的矩阵乘以$H$中系数的列向量之和为零的方程，再使用SVD算法找到$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]

H_from_points的函数首先检查输入点集的形状是否匹配，之后对输入点进行条件化。这是为了在数值计算过程中提高稳定性，尤其是在进行SVD分解时。条件化过程包括对每个点集计算平均值和标准差，然后用这些值创建一个条件矩阵。接下来，将条件矩阵与每个点集相乘，得到条件化后的点集。使用条件化后的点集创建线性方法所需的矩阵A。矩阵A包含两倍于对应点对数量的行，每行对应一个对应点对。对于每个对应点对，矩阵A的行存储了用于计算单应性矩阵H的信息。对矩阵A进行奇异值分解（SVD），得到U、S和V矩阵。从V矩阵中选择第8列（与单应性矩阵H的大小相匹配），并将其重塑为一个3x3矩阵。对H进行去条件处理，最后将归一化后的单应性矩阵H返回

3.1.2 仿射变换

仿射变换的6个自由度需要三个对应点来估计矩阵$H$，可以使用上述的DLT算法计算，也可以通过下述代码来计算

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]

Haffine_from_points的函数与前个函数前几步骤相同，在得到U、S和V矩阵后根据V矩阵的前两列创建B和C矩阵。根据B和C矩阵构建仿射变换矩阵H。再对H进行去条件处理，即将与条件化矩阵C1和C2相关的操作进行逆操作，最后将归一化后的仿射变换矩阵H返回

3.2 图像扭曲

对图像进行仿射变换就称为图像扭曲，一种简便的方法是使用Scipy中的ndimage包完成transformed_im = ndimage.affine_transform(im,A,b,size)线性变换A和平移向量b指定大小size对图像im应用仿射变换，主要代码如下

im = array(Image.open('filelist/PIL5.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]))


subplot(1, 2,  1)
gray()
imshow(im)
axis('off')

subplot(1, 2, 2)
gray()
imshow(im2)
axis('off')

show()

从效果上来看输出的图像结果中丢失的像素用零来填充

3.2.1 图像中的图像

将一副图像放在另一副图像中的某个位置就是仿射扭曲的一个例子，如下述代码，简单的把i m1放在im2的角点坐标上，先确定要填充的坐标在进行仿射变换

def image_in_image(im1,im2,tp):
    """ Put im1 in im2 with an affine transformation
        such that corners are as close to tp as possible.
        tp are homogeneous and counter-clockwise from top left. """ 
    
    # points to warp from
    m,n = im1.shape[:2]
    fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])
    
    # compute affine transform and apply
    H = homography.Haffine_from_points(tp,fp)
    im1_t = ndimage.affine_transform(im1,H[:2,:2],
                    (H[0,2],H[1,2]),im2.shape[:2])
    alpha = (im1_t > 0)
    
    return (1-alpha)*im2 + alpha*im1_t

在实际运用中，对于一般透视效果不明显的物体来说，一个仿射变换就足矣，但是在透视效果明显的情况下，同一个仿射变换就力不从心了。对于三个对应点，仿射变换可以进行图像扭曲，对应于六个自由度，因此可以将图像分成两个三角形分别进行图像扭曲使得在透视明显的情况下也能达到较好的效果，实现代码和效果如下

im2 = array(Image.open('filelist/PIL5.jpg').convert('L'))
im1 = array(Image.open('filelist/tp1.jpg').convert('L'))

m, n = im1.shape[:2]
fp = array([[0, m, m, 0], [0, 0, n, n], [1, 1, 1, 1]])
tp = array([[1895,2208,2260,1881],[1213,1151,674,694],[1,1,1,1]])

tp2 = tp[:,:3]
fp2 = fp[:,:3]

# compute affine transform and apply
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 = warp.alpha_for_triangle(tp2,im2.shape[0],im2.shape[1])
im3 = (1-alpha)*im2 + alpha*im1_t

tp2 = tp[:,[0,2,3]]
fp2 = fp[:,[0,2,3]]

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 = warp.alpha_for_triangle(tp2,im2.shape[0],im2.shape[1])
im4 = (1-alpha)*im3 + alpha*im1_t

gray()
imshow(im4)
axis('off')

show()

和之前的图像相比，这里插入的图像边缘更加平滑，有透视的效果

3.2.2 分段放射扭曲

三角形图像块的仿射扭曲可以完成角点的精确匹配，对应点对集合间常用的扭曲方式为分段仿射扭曲。给定图像的标记点，对它们进行三角剖分并使用仿射扭曲扭曲每个三角形，就能实现两幅图像对应点的扭曲对应。三角化点通常使用洛克三角剖分法

from scipy.spatial import Delaunay

x,y = array(random.standard_normal((2,100)))
tri = Delaunay(np.c_[x, y]).simplices

figure()
for t in tri:
    t_ext = [t[0], t[1], t[2], t[0]] # add first point to end
    plot(x[t_ext],y[t_ext],'r')

plot(x,y,'*')
axis('off')
show()

结果如上，随机生成数据实例点，每三个点构成一个三角形，且三角剖分中所有三角形的最小角度最大，所有点连接起来形成一个网络。在实际运用中，就是将图像中的某些目标点选取出来，将一副图像扭曲到另一副图像的非平坦区域

3.2.3 图像配准

图像配准主要目标是将多张不同时间、不同视角、不同传感器拍摄的图像进行对齐，使得它们在空间上能够完美地对应起来。它的主要步骤如下：

1. 特征提取：从输入图像中提取出可以表示图像特征的关键点或区域。常用的特征提取方法有 SIFT、SURF、ORB、FAST 等
2. 特征匹配：根据提取的特征，在两幅图像之间找到特征之间的对应关系。这一步通常使用最近邻搜索或基于特征描述子的匹配方法
3. 建立变换模型：根据匹配的特征点对，建立一个描述图像间变换关系的模型。常用的变换模型有平移、旋转、仿射变换等，也可以使用更复杂的变换模型，如非线性变换或仿射加旋转的复合变换
4. 优化变换参数：求解变换模型的参数，使得两张图像间的误差（如欧几里得距离、互信息等）最小化
5. 插值和融合：根据求解出的变换参数，对输入图像进行插值和融合，得到配准后的图像

书上的实验以一个人在一年中每天的人脸进行图像配准，其实就是将人脸的各个特征部位对齐，使人脸的尺寸、方向和位置都统一起来，由于没有人脸的图像文件，此处省略实验

3.3 创建全景图

创建全景图像是指将一系列具有重叠区域的局部图像拼接成一张完整的全景图像，在我们手机的照相机里上都有这个功能，更为优秀的是，除了能达到平面的图像展示外还能以3D的全景图像方式展示，它包含以下基本步骤：

1. 采集图像：使用具有固定焦距的相机拍摄一系列局部图像
2. 特征提取与匹配：从局部图像中提取特征点（如 SIFT、SURF、ORB 等），并在不同局部图像之间找到特征之间的对应关系
3. 姿态估计：根据匹配的特征点对，计算局部图像之间的相对旋转和平移关系
4. 全景图像拼接：将局部图像根据估计的变换关系进行拼接
5. 融合与优化：对拼接后的图像进行融合和优化处理，消除拼接缝、校正畸变等

3.3.1 RANSAC

RANSAC（Random Sample Consensus，随机抽样一致） 主要目标是从一组包含异常值的数据中，找到一个合适的模型来描述数据的内在规律。RANSAC 算法的核心思想是通过随机抽样的方式，迭代地寻找数据中的“局内点”（Inliers，即符合模型的数据点），并使用这些局内点来计算最优模型参数

RANSAC 算法的优点是能够有效地处理含有噪声的数据，并且对于局内点的比例要求不高。然而，RANSAC 算法也有一些缺点，如计算复杂度较高、对局内点和局外点的比例敏感等

以下是一个简单的 Python 示例，使用 RANSAC 算法进行直线拟合。在这个例子中，我们假设有一组包含噪声的二维数据点，我们希望找到一个最优直线来描述这些数据点。

3.3.2 拼接图像

创建全景图像的步骤在之前已经提过，已经解决之前的.sift文件生成的问题，现在使用以下代码进行全景拼接

featname = ['./filelist/tian' + str(i + 1) + '.sift' for i in range(3)]
imname = ['./filelist/tian' + str(i + 1) + '.jpeg' for i in range(3)]

# extract features and match
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])


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


# 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)

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

该图像一共由三张图像拼接而成，但是从结果来看，左侧和中间的图像拼接的并不好，可能是由于两张图像的特征匹配不足。中间和右侧的图片匹配程度较好，但还是能看到有岷县的分割线