计算机视觉课后作业3——图像变换

                               图像变换

一.实验内容：

1. H_from_points(fp,tp)和affine_from_points函数，根据匹配点，通过最小二乘法分别求解单应性变换矩阵和仿射变换矩阵
2. 根据仿射或单应性变换实现图像的扭曲、映射、融合

二.实验平台：PyCharm

三.具体实验：

1. 单应性变换——H_from_points(fp,tp)函数和affine_from_points函数

本实验中单应性变换主要涉及直接线性变换算法和仿射变换

简介：

单应性变换其实就是一个平面到另一个平面的变换关系，是将一个平面内的点映射到另一个平面内的二维投影变换。单应性变换通过找到图片中相同的物理点，也就是对应点，将一个图像中的点映射到另一个图像中的对应点。单应性变换具有很强的实用性，比如图像配准，图像纠正和纹理扭曲，以及创建全景图像，都会使用单应性变换。

实验原理：

a.首先我们需要了解齐次坐标的概念。所谓齐次坐标就是将一个原本是n维的向量用一个n+1维向量来表示。例如，二维点(x,y)的齐次坐标表示为(hx,hy,h)。点的齐次坐标是依赖于其尺度定义的，所以，x=[x,y,w]=[ax,ay,aw]=[x/w,y/w,1]都表示同一个二维点。

b.单应性变换H，按照下面的方程映射二维中的点（齐次坐标意义下）：

 

在a中提到的，点的点的齐次坐标是依赖于其尺度定义的。因此，单应性矩阵H也仅依赖尺度定义，所以，单应性矩阵具有8个独立的自由度。我们通常使用w=1来归一化点，这样，点具有唯一的图像坐标x和y。这个额外的坐标是的我们可以简单地使用一个矩阵来表示变换。

c.下面的函数可以实现对点进行归一化和转换齐次坐标的功能：

def normallize(points):
    """在齐次坐标意义下，对点集进行归一化，是最后一行为1"""
    for row in points:
        row /= points[-1]
    return points

def make_homog(points):
    """将点集（dim×n的数组）转换为齐次坐标表示"""

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

进行点和变换的处理时，我们会按照列优先的原则存储这些点。因此，n个二维点集将会存储为齐次坐标意义下的一个3×n数组。

直接线性变换算法：

原理：

单应性矩阵可以通过两幅图像（或者平面）中对应点对计算出来。一个完全射影变换具有8个自由度。根据对应点约束，每个对应点对可以写出两个方程，分别对应于x和y坐标。因此，计算单应性矩阵H需要4个对应点对。

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

函数代码如下：

def H_from_points(fp, tp):
    """使用线性DLT方法，计算单应性矩阵H，使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]

仿射变换：

原理：

是一种二维坐标到二维坐标之间的线性变换，它保持了二维图形的“平直性”（直线经过变换之后依然是直线）和“平行性”（二维图形之间的相对位置关系保持不变，平行线依然是平行线，且直线上点的位置顺序不变），但是角度会改变。任意的仿射变换都能表示为乘以一个矩阵(线性变换)，再加上一个向量 (平移) 的形式。 

仿射变换具有6个自由度，因此我们需要三个对应点来估计矩阵H。通过将最后两个元素设置为0，即h7=h8=0。

函数代码如下：

    """计算H仿射变换，使得tp是fp经过仿射变换H得到的"""

    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() # 两个点集，必须都进行相同的缩放
    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)
    
    # 如Hartley和Zisserman著的Multiplr View Geometry In Computer,Scond Edition所示，
    # 创建矩阵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]　

2.图像扭曲：

对图像块应用仿射变换，我们将其称为图像扭曲（或者仿射扭曲）。该操作不仅经常在计算机图形学中，而且经常出现在计算机视觉算法总。扭曲的操作可以使用SciPy工具包中的ndimage包来简单完成。

代码：

from matplotlib.pyplot import *
from scipy import ndimage
from PIL import Image
im = array(Image.open(r'C:\Users\oulia\Pictures\b1.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]))

gray()
subplot(121)
imshow(im)
axis('off')
subplot(122)
imshow(im2)
axis('off')
show()

实验结果截图：

3.图像融合：

ps:因为书本上的代码需要修改比较多的库，所以没有采用，使用的是PIL库中的paste方法。

步骤：
1）从图像中复制子矩形

2）处理一个子矩形，并将其粘贴回去

注意：当粘贴区域返回时，区域的大小必须与给定区域完全匹配。此外，该区域不能扩展到图像之外。

代码：

from PIL import Image
base_img = Image.open(r'C:\Users\oulia\Pictures\c3.jpg')
tmp_img = Image.open(r'C:\Users\oulia\Pictures\a1.jpg')
box = (0, 0, 600, 800)  # 要替换的区域
region = tmp_img  # 整图替换
region = region.resize((box[2] - box[0], box[3] - box[1]))
base_img.paste(region, box)
base_img.show()

实验结果截图：