文章目錄

• 圖像到圖像的映射
• • 1 單應性變換
  • • 1.1 直接線性變換算法
    • 1.2 仿射變換
  • 2 圖像扭曲
  • • 2.1 圖像中的圖像
    • 2.2 分段仿射扭曲
    • 2.3 圖像配準
  • 3 創(chuàng)建全景圖
  • • 3.1 RANSAC
    • 3.2 穩(wěn)健的單應性矩陣估計
    • 3.3 拼接圖像

圖像到圖像的映射

1 單應性變換

單應性變換是將一個平面內的點映射到另一個平面內的二維投影變換。

單應性變換具有很強的實用性，比如圖像配準、圖像糾正和紋理扭曲，以及創(chuàng)建全景圖像。

本質上，單應性變換H，按照下邊的方程映射二維中的點（齊次坐標意義下）：

對于圖像平面內的點，齊次坐標是一個非常有用的表示方式。點的齊次坐標是依賴于其尺度定義的，所以x=[x, y, w]=[αx, αy, αw]=[x/w, y/w,1]都表示同一個二維點。因此，單應性矩陣H也僅僅依賴尺度定義，所以單應性矩陣具有9個獨立的自由度。

通常使用w=1來歸一化，這樣點具有唯一的圖像坐標x和y，這樣可以簡單地使用一個矩陣來表示變換。

創(chuàng)建homography.py文件，使用下邊的函數(shù)實現(xiàn)對點進行歸一化和轉換齊次坐標的功能：

def normalize(points):for row in points:row /= points[-1]return pointsdef make_homog(points):return vastack((points,ones((1,points,shape[1]))))

行點和變換處理時，我們會按照列優(yōu)先的原則存儲這些點。因此，n個二維點集將會存儲為齊次坐標意義下的一個3 *n數(shù)組。這種格式使得矩陣乘法和點的變換操作更加容易。

在這些投影的變換中，有一些特別重要的變換。比如，仿射變換可以應用于圖像扭曲、相似變換可以應用于圖像配準。

1.1 直接線性變換算法

單應性矩陣可以由兩幅圖像中對應點對計算出來。一個完全射影變換具有8個自由度，根據對應點約束，每個對應點對可以寫出兩個方程，分別對應于x和y坐標。因此，計算單應性矩陣H需要4個對應點對。

DLT（直接線性變換）是給定4個或者更多對應點對矩陣，來計算單應性矩陣H的算法。將單應性矩陣H作用于對應點對上，重新寫出一個齊次方程Ah=0，其中A是一個具有對應點對二倍數(shù)量行數(shù)的矩陣。將這些對應點對方程的系數(shù)堆疊到一個矩陣中，可以使用奇異值分解找到H的最小二乘解。

下面是實現(xiàn)該算法的代碼：

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-9C1 = diag([1/maxstd,1/maxstd,1])C1[0][2] = -m[0]/maxstdC1[1][2] = m[1]/maxstdfp = dot(C1,fp)m = mean(tp[:2],axis=1)maxstd = max(std(tp[:2],axis=1)) +1e-9C2 = diag([1/maxstd,1/maxstd,1])C2[0][2] = -m[0]/maxstdC2[1][2] = -m[1]/maxstdtp = 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[0][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]

上面函數(shù)的第一步操作是檢查點對的兩個數(shù)組中點的數(shù)目是否相同。如果不相同，函數(shù)拋出異常信息。

對這些點進行歸一化操作，使其均值為0，方差為1。然后使用對應點對構造矩陣A。最小二乘解即為矩陣SVD分解后所得矩陣V的最后一行。該行經過變形后得到矩陣H。然后對著矩陣進行處理和歸一化，返回輸出。

1.2 仿射變換

由于仿射變換具有6個自由度，因此我們需要三個對應點對來估計矩陣H。通過將最后兩個元素設置為0，即h7 = h8 = 0。仿射變換可以用上面的DLT算法估計得出。

下面的函數(shù)使用對應點對來計算仿射變換矩陣：

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-9C1 = diag([1/maxstd,1/maxstd,1])C1[0][2] = -m[0]/maxstdC1[1][2] = -m[1]/maxstdfp_cond = dot(C1,fp)m = mean(tp[:2],axis=1)C2 = C1.copy() C2[0][2] = -m[0]/maxstdC2[1][2] = -m[1]/maxstdtp_cond = dot(C2,fp)A = concatenate((fp_cond[:2],tp_cond[:2]),axis=0)U,S,V = linalg.svd(A.T)tmp = V[:2].TB = 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包來簡單完成。命令：

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

使用一個線性變換A和一個平移向量b來對圖像塊應用仿射變換。size用來指定輸出圖像的大小。默認輸出圖像設置為和原始圖像同樣大小。

了研究該函數(shù)是如何工作的，我們可以運行下面的命令：

from PIL import Image from pylab import * from scipy import ndimageim = array(Image.open('jimei2.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(121) gray() axis('off') imshow(im)subplot(122) gray() axis('off') imshow(im2)show()

得到的結果如下圖所示，可見輸出圖像中丟失的像素用零來填充：

圖1 用仿射變換扭曲圖像

2.1 圖像中的圖像

仿射扭曲的一個簡單例子是將圖像或者圖像的一部分放置在另一幅圖像中，使得他們能夠和指定的區(qū)域或者標記物對齊。

以下函數(shù)的輸入參數(shù)為兩幅圖像和一個坐標。該坐標為將一副圖像放置到第二幅圖象中的角點位置：

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 = Haffine_form_points(tp,fp)im1_t = ndimage.affine_transform(im,H[:2,:2],(H[0,2],H[1,2]),im2.shape[:2])alpha = (im1_t>0)return (1-alpha)*im2 + alpha*im1_t

將扭曲的圖像和第二幅圖像融合，就創(chuàng)建alpha圖像。該圖像定義了每個像素從各個圖像中獲取的像素值成分多少。這里基于以下事實：扭曲的圖像是在扭曲區(qū)域邊界之外以0來填充的圖像，來創(chuàng)建一個二值的alpha圖像。嚴格意義上，需要在第一幅圖象中的潛在0像素上加上一個小的數(shù)值，或者合理的處理這些0像素。

下面幾行代碼會將一幅圖像插入另一幅圖像：

im4= array(Image.open('jimei.jpg').convert('L')) im5 = array(Image.open('guanggaopai.jpg').convert('L'))tp = array([[50,300,300,50],[10,10,600,550],[1,1,1,1]]) im6 = image_in_image(im4,im5,tp) figure() gray() imshow(im6) axis('equal') axis('off') show()

圖2 用仿射變換將一幅圖像放置到另一幅圖像中

函數(shù)Haffine_from points() 會返回給定對應點對的最有仿射變換。在上面的例子中，對應點對為圖像和公告牌的角點。如果透視效應比較弱，那么這種方法會返回很好的結果。

2.2 分段仿射扭曲

給定任意圖像的標記點，通過將這些點進行三角剖分，然后使用仿射扭曲來扭曲每個三角形，我們可以將圖像和另一幅圖像的對應標記點扭曲對應。為了三角化這些點，我們經常使用狄洛克三角剖分方法。在matplotlib中有狄洛克三角剖分：

from PIL import Image from pylab import * import numpy as np from scipy.spatial import Delaunayx,y = array(np.random.standard_normal((2,100))) """centers, edges, tri, neighbors = md.delaunay(x, y""" tri = Delaunay(np.c_[x,y]).simplicesfigure()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()

圖3 隨機二維點集的狄洛克三角
現(xiàn)在讓我們將該算法用于一個例子，在該例子中，在5*6的網格上使用30個控制點，將一幅圖像扭曲到另一幅圖像中的非平坦區(qū)域。

首先，我們需要寫出一個用于分段仿射圖像扭曲通用的扭曲函數(shù)。下面的代碼實現(xiàn)該功能：

def pw_affine(fromim,toim,fp,tp,tri):im = toim_copy()is_color = len(fromim.shape) == 3 im_t = zeros(im.shape,'unit8')for t in tri:H = Haffine_from_points(tp[:,t],fp[:,t])if is_color:for col in range(fromim,shape[2]):im_t[:,:,col] = ndimage.affine_transform(fromim[:,:,col],H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])else:im_t = ndimage.affine_transform(fromim,H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])alpha = alpha_from_triangle(tp[:,t],im.shape[0],im.shape[1])im[alpha>0] = im_t[alpha>0]return im

在該代碼中，我們首先見檢查該圖像是灰度圖像還是彩色圖像。如果圖像為彩色圖像，則對每個顏色通道進行扭曲處理。因為對于每個三角形來說，仿射變換是唯一確定的，所以我們這里使用Haffine_from_points()函數(shù)來處理。

2.3 圖像配準

圖像配準是對圖像進行變換，使得變換后的圖像能夠在常見的坐標系中對齊。配準可以是嚴格配準，也可以是非嚴格配準。為了能夠進行圖像對比和更精細的圖像分析，圖像配準是一步非常重要的操作。

目前，很難找到一種普適的方法能夠應對所有的配準情況，任何一種配準算法都必須考慮圖像的成像原理、幾何變形、噪聲影響、配準精度等因素。不過，從原理上將，配準算法可以大致分為以下四個步驟：

（1）特征提取
采用人工或者自動的方法檢測圖像中的不變特征，如：閉合區(qū)域、邊緣、輪廓、角點等。特征提取算法需要滿足三個條件
（a）顯著性，所提取的特征應該是比較明顯的，分布廣泛的、易于提取的特征；
（b）抗噪性，具有較強的噪聲抑制能力且對成像條件的變化不敏感；
（c）一致性，能準確地檢測出兩幅圖像的共有特征；

（2）特征匹配
通過特征描述算作及相似性度量來建立所提取的特征之間的對應關系。特征匹配常用到的區(qū)域灰度、特征向量空間分布和特征符號描述等信息。某些算法在進行特征匹配的同時也完成了變換模型參數(shù)的估計；

（3）變換模型估計
指根據待配準圖像與參考圖像之間的幾何畸變的情況，選擇能最佳擬合兩幅圖像之間變化的幾何變換模型，可以分為全局映射模型和局部映射模型。其中，全局映射模型利用所有控制點信息進行全局參數(shù)估計；局部映射模型利用圖像局部的特征分別進行局部參數(shù)估計。常見的變換模型包括仿射變換、透視變換、多項式變換等，其中最常用的是仿射變換和多項式變換。

（4）坐標變換與插值
將輸入圖像做對應的參數(shù)變換，使它與參考圖像處于同一個坐標系下。由于圖像變換后的坐標點不一定是整數(shù)，因此，需要考慮一定的插值處理操作。常用的插值方法包括：最近鄰插值、雙線性插值、雙三次插值、B樣條插值、高斯插值。

3 創(chuàng)建全景圖

3.1 RANSAC

RANSAC是用來找到正確模型來擬合帶有噪聲數(shù)據的迭代方法。給定一個模型，例如點集之間的單應性矩陣，該算法的基本思想是數(shù)據中包含正確的點和噪聲點，合理的模型能夠在描述正確點的同時摒棄噪聲點。

3.2 穩(wěn)健的單應性矩陣估計

在使用RANSAC模塊時，需要在相應的類中實現(xiàn)fit()和get_error()方法，剩下的就是正確使用ransac.py。

將下面的代碼添加到homography.py中：

class RansacModel(object): def __init__(self,debug=False):self.debug = debugdef fit(self, data):data = data.Tfp = data[:3,:4]tp = data[3:,:4]return H_from_points(fp,tp)def get_error( self, data, H):data = data.Tfp = data[:3]tp = data[3:]fp_transformed = dot(H,fp)fp_transformed = normalize(fp_transformed)return sqrt( sum((tp-fp_transformed)**2,axis=0) )

可以看到，這個類包含fit() 方法。該方法僅僅接受由ransac.py選擇的4個對應點對，然后你和一個單應性矩陣。記住，4個點對是計算單應性矩陣所需的最少數(shù)目。由于get_error()方法對于每個對應點對使用該單應性矩陣，然后返回相應的平方距離之和。因此ransac算法能夠判定哪些點對是正確的，那些是錯誤的。在實際中，我們需要在距離上使用一個闕值來決定哪些矩陣是合理的。為了使用方便，將下面的函數(shù)添加到homography.py 文件中：

def H_from_ransac(fp,tp,model,maxiter=1000,match_theshold=10):import ransacdata = vstack((fp,tp))H,ransac_data = ransac.ransac(data.T,model,4,maxiter,match_theshold,10,return_all=True)return H,ransac_data['inliers']

該函數(shù)同樣允許提供闕值和最小期望的點對數(shù)目。最重要的參數(shù)是最大迭代次數(shù)：程序推出太早可能得到一個壞解；迭代次數(shù)太多會占用太多空間。函數(shù)的返回結果為單應性矩陣和對應該單應性矩陣的正確點對。類似于下面的操作，你可以將RANSAC算法應用于點對上：

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,tpmodel = 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)[0]

在該例子中，圖像2是中心圖像，也是我們希望將其他圖像變成的圖像。圖像0和圖像1應該是從右邊扭曲，圖像3和圖像4應該是從左邊扭曲。在每個圖像對中，由于匹配是從最右邊的圖像計算出來的。所以我們將對應的順序進行了顛倒，使其從左邊圖像開始扭曲。因為我們不關心該扭曲例子中的正確點對，所以僅需要該函數(shù)的第一個輸出。

3.3 拼接圖像

估計出圖像間的單應性矩陣，現(xiàn)在我們需要將所有圖像扭曲到一個公共的圖像平面上。通常，這里的公共平面為中心圖像平面。一種方法是創(chuàng)建一個很大的圖像，比如圖像中全部填充0，使其和中心圖像平行，然后將所有的圖像扭曲到上面。由于我們所有的圖像是由照相機水平旋轉拍攝的，因此我們需要一個比較簡單的步驟：將中心圖像左邊或者右邊的區(qū)域填充0，以便為扭曲的圖像騰出空間。

def panorama(H,fromim,toim,padding=2400,delta=2400):is_color = len(fromim.shape) == 3def transf(p):p2 = dot(H,[p[0],p[1],1])return (p2[0]/p2[2],p2[1]/p2[2])if H[1,2]<0: print('warp - right')if is_color: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: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)if is_color: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:toim_t = hstack((zeros((toim.shape[0],padding)),toim))fromim_t = ndimage.geometric_transform(fromim,transf,(toim.shape[0],toim.shape[1]+padding))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

對于通用的geometric_transform()函數(shù)，我們需要指定能夠描述像素到像素間映射的函數(shù)。在這個例子中，transf()函數(shù)就是該指定函數(shù)。干函數(shù)通過將像素和H相乘，然后對齊次坐標進行歸一化來實現(xiàn)像素間的映射。通過查看H中的平移量，我們可以決定應該將圖像填補到左邊還是右邊。當圖像填補到左邊時，由于目標圖像中點的坐標也變化了，所以在“左邊“情況中，需要在單應性矩陣中加入平移。簡單起見，我們同樣使用0像素的技巧來尋找alpha圖。

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)figure() imshow(array(im_42)) axis('off') show()

圖4使用SIFT對應點自動創(chuàng)建水平全景圖像

總結

以上是生活随笔為你收集整理的计算机视觉编程——图像到图像的映射的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。