最近一直在學習Python計算機視覺程式設計中影象到影象之間的對映，這些變化可以用於影象扭曲變形和影象配準。所謂單應性變換就是將一個平面內的點對映到另一個平面內的二維投影變換。單應性變換具有很強的實用性，比如影象配準、影象糾正和紋理扭曲，以及建立全景影象等。其實，單應性變換H，是按照下面的公式對映二維中的點（其次座標意義下）：

齊次座標是在實際應用中非常有用的表示式，點的齊次座標依賴其尺度定義，所以x=[x,y,w]=[ax,ay,aw]=[x/w,y/w,1]都可以表示同一個二維點，所以單應性矩陣H也僅僅依賴尺度定義，單應性矩陣H有8個獨立的自由度。我們通常使用w=1來歸一化點，這樣點具有唯一的影象座標x和y。

下面的程式碼實現了對點進行歸一化和轉換成齊次座標：

# coding=UTF-8
from PIL import Image
from numpy import *
from pylab import *


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


def normalize(points):
    '''
    在齊次座標下進行歸一化處理
    :param points:
    :return:
    '''
    for row in points:
        row /= points[-1]
    return points
 

def Hfrompoints(fp, tp):
    '''
    求解單應性矩陣，使fp對映到tp
    :param fp:
    :param tp:
    :return:
    '''
    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[0] / 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[0] / maxstd
    tp = dot(C1, tp)

    # 求解 H
    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]