计算机视觉（三）：图像到图像的映射

这一次主要是讲图像之间的变换，以及一些计算变换的实用方法。这些变换可以用于图像扭曲变形和图像配准。

一、单应性变换

单应性变换： 将一个平面内的点映射到另一个平面内的二维投影变换。平面指图像或者三维中的平面表面。
单应性变换H按照下面的方程映射二维中的点：

齐次坐标： 点的齐次坐标依赖于其尺度定义的所以以下公式都表示同一个二维点

单映性矩阵H也依赖尺度定义，单应射矩阵具有8个独立的自由度，使用w=1来归一化点，这样，点具有唯一的图像坐标x和y。
我们可以使用以下函数实现归一化和转换齐次坐标的功能。在进行点和变换的处理时，我们会按照列优先的原则存储这些点。因此，n个二维点集将会存储为齐次坐标意义下的一个3xn数组

from numpy import *

def normalize(points):
    for row in points:
        row /= points[-1]
    return points


def make_homog(points):
    return vstack((points ,ones((1 ,points.shape[1]))))

1.1 直线线性变换算法

单应性矩阵可以有两幅图像中对应点对计算出来，一个完全射影变换具有8个自由度，所以计算一个单应性矩阵H需要4个对应点对。
DLT（直接线性变换） 是给定4个或者更多对应点对矩阵，来计算单应性矩阵H的算法。将单应性矩阵H作用在对应点对上，我们可以得：

或者Ah=0，A是一个具有对应点对二倍数量行数得矩阵。将这些对应点对方程的系数对堆叠到一个矩阵中，可以使用SVD(奇异值分解)算法找到H的最小二乘解。代码如下：

def H_from_points(fp,tp):
    if fp.shape != tp.shape:
        raise RuntimeError('number of points do not match')
        
    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)

    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)
    
    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))    
    
    H = dot(linalg.inv(C2),dot(H,C1))
    
    return H / H[2,2]

对点进行归一化操作，使其均值为0，方差为1。算法的稳定性取决于坐标的表示情况和部分数值计算的问题，所以归一化操作非常重要！！！

1.2 仿射变换

仿射变换具有六个自由度，所以我们需要三个对应点对来估计矩阵H。通过将最后两个元素设置为0，即7h7=h8=0，仿射变换可以用上面的DLT算法估计得出。

但是我们还有另外一种计算单应性矩阵H的算法，如下所示：

def Haffine_from_points(fp,tp):
    
    if fp.shape != tp.shape:
        raise RuntimeError('number of points do not match')
        
    # 映射起始点
    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)
    
    # 映射对应点
    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)
    
    # 归一化后点的均值为0，所以平移量为0
    A = concatenate((fp_cond[:2],tp_cond[:2]), axis=0)
    U,S,V = linalg.svd(A.T)
    
    # 创建B.C矩阵
    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]))
    
    # 反归一化
    H = dot(linalg.inv(C2),dot(H,C1))
    
    return H / H[2,2]

二、图像扭曲

对图像块应用仿射变换也就是图像扭曲或者仿射扭曲。
扭曲操作命令如下：

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

演示代码及演示结果如下所示：

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

im = array(Image.open('yellow.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()
imshow(im2)
show()

效果图：

原图：

2.1 图像中的图像

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

可用如下函数实现，输入参数为两幅图像和一个坐标，坐标为第一幅图像放置到第二幅图像中的交点坐标：

def image_in_image(im1 ,im2 ,tp):

    m ,n = im1.shape[:2]
    fp = array([[0 ,m ,m ,0] ,[0 ,0 ,n ,n] ,[1 ,1 ,1 ,1]])

    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

创建alpha图像，该图像定义了每个像素从各个图像中获取的像素值成分多少，扭曲的图像实在扭曲区域边界之外以0来填充的图像，来创建一个alpah图像。

仿射变换实现：

import warp
im1 = array(Image.open('yellow.jpg').convert('L'))
im2 = array(Image.open('1.jpg').convert('L'))
tp = array([[264,538,540,264],[40,36,605,605],[1,1,1,1]])
im3 = warp.image_in_image(im1,im2,tp)
figure()
gray()
imshow(im3)
axis('equal')
axis('off')
show()

2.2 分段仿射扭曲

给定任意图像的标记点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，可以将图像和另一幅图像的对应标记点扭曲对应。
为了三角化这些点，我们经常使用狄洛克三角剖分 方法。在Matplotlib （但是不在 PyLab 库中）中有狄洛克三角剖分，可以用下面的方式使用它：

import matplotlib.delaunay as md
x,y = array(random.standard_normal((2,100)))
centers,edges,tri,neighbors = md.delaunay(x,y)
figure()
for t in tri:
 t_ext = [t[0], t[1], t[2], t[0]] # 将第一个点加入到最后
 plot(x[t_ext],y[t_ext],'r')
plot(x,y,'*')
axis('off')
show()

三角剖分函数：

import matplotlib.delaunay as md
def triangulate_points(x,y):
 """ 二维点的Delaunay 三角剖分"""
 centers,edges,tri,neighbors = md.delaunay(x,y)
 return tri

2.3 图像配准

图像配准 对图像进行变换，使变换后的图像能够在常见的坐标系中对齐。配准可以是严格配准，也可以是非严格配准。

2.3.1 特征提取

特征提取是指分别提取两幅图像中共有的图像特征，这种特征是出现在两幅图像中对比列、旋转、平移等变换保持一致性的特征，如线交叉点、物体边缘角点、虚圆闭区域的中心等可提取的特征。特征包括：点、线和面三类

三、创建全景图

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

3.1 RANSAC

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

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

局外点产生的原因：噪声的极值；错误的测量方法；对数据的错误假设。

3.2 稳健的单应性矩阵估计

在使用 RANSAC 模块时，我们只需要在相应 Python 类中实现 fit() 和 get_error()方法，剩下就是正确地使用 ransac.py。
使用SIFT特征自动找到匹配对应，调用函数如下：

**import sift
featname = ['Univ'+str(i+1)+'.sift' for i in range(5)]
imname = ['Univ'+str(i+1)+'.jpg' for i in range(5)]
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])**

3.3 拼接图像

估计出图像间的单应性矩阵（使用 RANSAC 算法），现在需要将所有的图像扭曲到一个公共的图像平面上，可以使用一个较简单的步骤：将中心图像左边或者右边的区域填充 0，以便为扭曲的图像腾出空间：

def panorama(H,fromim,toim,padding=2400,delta=2400):
 # 检查图像是灰度图像，还是彩色图像
 is_color = len(fromim.shape) == 3
 # 用于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 在右边
 print 'warp - right'
 # 变换fromim
 if is_color:
 # 在目标图像的右边填充0
 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],
 transf,(toim.shape[0],toim.shape[1]+padding))
 else:
 # 在目标图像的右边填充0
 toim_t = hstack((toim,zeros((toim.shape[0],padding))))
 fromim_t = ndimage.geometric_transform(fromim,transf,
 (toim.shape[0],toim.shape[1]+padding))
 else:
 print 'warp - left'
 # 为了补偿填充效果，在左边加入平移量
 H_delta = array([[1,0,0],[0,1,-delta],[0,0,1]])
 H = dot(H,H_delta)
 # fromim 变换
 if is_color:
 # 在目标图像的左边填充0
 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],
 transf,(toim.shape[0],toim.shape[1]+padding))
 else:
 # 在目标图像的左边填充0
 toim_t = hstack((zeros((toim.shape[0],padding)),toim))
 fromim_t = ndimage.geometric_transform(fromim,
 transf,(toim.shape[0],toim.shape[1]+padding))
 # 协调后返回（将fromim 放置在toim 上）
 if is_color:
 # 所有非黑色像素
 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)
 return toim_t

四、校园全景拼接图

from cmath import sqrt

import numpy as np
from PCV.geometry import warp
from PIL.Image import Image
from numpy import dot, hstack, zeros, array, vstack, argsort, arccos
from scipy import ndimage, linalg

from code.second import imname

from code.homography import normalize

from code import homography


class RansacModel(object):
    """ Class for testing homography fit with ransac.py from
        http://www.scipy.org/Cookbook/RANSAC"""

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

    def fit(self, data):
        """ Fit homography to four selected correspondences. """

        # transpose to fit H_from_points()
        data = data.T

        # from points
        fp = data[:3, :4]
        # target points
        tp = data[3:, :4]

        # fit homography and return
        return H_from_points(fp, tp)

    def get_error(self, data, H):
        """ Apply homography to all correspondences,
            return error for each transformed point. """

        data = data.T

        # from points
        fp = data[:3]
        # target points
        tp = data[3:]

        # transform fp
        fp_transformed = dot(H, fp)

        # normalize hom. coordinates
        fp_transformed = normalize(fp_transformed)

        # return error per point
        return sqrt(sum((tp - fp_transformed) ** 2, axis=0))

def H_from_ransac(fp, tp, model, maxiter=1000, match_theshold=10):
    """ Robust estimation of homography H from point
        correspondences using RANSAC (ransac.py from
        http://www.scipy.org/Cookbook/RANSAC).

        input: fp,tp (3*n arrays) points in hom. coordinates. """

    import ransac

    # group corresponding points
    data = vstack((fp, tp))

    # compute H and return
    H, ransac_data = ransac.ransac(data.T, model, 4, maxiter, match_theshold, 10, return_all=True)
    return H, ransac_data['inliers']

def match(desc1, desc2):
    desc1 = array([d / linalg.norm(d) for d in desc1])
    desc2 = array([d / linalg.norm(d) for d in desc2])

    dist_ratio = 0.6
    desc1_size = desc1.shape

    matchscores = zeros((desc1_size[0]), 'int')
    desc2t = desc2.T  # 预先计算矩阵转置
    for i in range(desc1_size[0]):
        dotprods = dot(desc1[i, :], desc2t)  # 向量点乘
        dotprods = 0.9999 * dotprods
        # 反余弦和反排序，返回第二幅图像中特征的索引
        indx = argsort(arccos(dotprods))

        # 检查最近邻的角度是否小于 dist_ratio 乘以第二近邻的角度
        if arccos(dotprods)[indx[0]] < dist_ratio * arccos(dotprods)[indx[1]]:
            matchscores[i] = int(indx[0])

    return matchscores

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)
    return fp, tp


model = homography.RansacModel()

fp, tp = convert_points(1)
H_12 = homography.H_from_ransac(fp, tp, model)[0]

fp, tp = convert_points(0)
H_01 = homography.H_from_ransac(fp, tp, model)[0]

tp, fp = convert_points(2)
H_32 = homography.H_from_ransac(fp, tp, model)[0]

tp, fp = convert_points(3)
H_43 = homography.H_from_ransac(fp, tp, model, match_threshold=3000)[0]

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],
                                                                  transf, (toim.shape[0], toim.shape[1] + padding))
        else:
            # pad the destination image with zeros to the right
            toim_t = hstack((toim, zeros((toim.shape[0], padding))))
            fromim_t = ndimage.geometric_transform(fromim, transf,
                                                   (toim.shape[0], toim.shape[1] + padding))
    else:
        print('warp - left')
        # add translation to compensate for padding to the 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],
                                                                  transf, (toim.shape[0], toim.shape[1] + padding))
        else:
            # pad the destination image with zeros to the left
            toim_t = hstack((zeros((toim.shape[0], padding)), toim))
            fromim_t = ndimage.geometric_transform(fromim,
                                                   transf, (toim.shape[0], toim.shape[1] + padding))

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

    return toim_t


im1 = np.array(Image.open(imname[1]))
delta = im1.shape[1]  # 領域追加と水平移動量

im2 = np.array(Image.open(imname[2]))
im_12 = warp.panorama(H_12, im1, im2, delta, delta)

im1 = np.array(Image.open(imname[0]))
im_02 = warp.panorama(np.dot(H_12, H_01), im1, im_12, delta, delta)

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

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