【python计算机视觉编程——3.图像到图像的映射】

python计算机视觉编程——3.图像到图像的映射

3.图像到图像的映射
- 3.1 单应性变换
- - 3.1.1 直接线性变换算法（DLT）
  - 3.1.2 仿射变换
- 3.2 图像扭曲
- - 3.2.1 图像中的图像
  - 3.2.2 分段仿射扭曲
  - 3.2.3 图像配准
- 3.3 创建全景图
- - 3.3.1 RANSAC
  - 3.3.2 稳健的单应性矩阵估计
  - 3.3.3 拼接图像

3.图像到图像的映射

3.1 单应性变换

将一个图像平面上的点映射到另一个图像平面上的技术。它在图像配准、立体视觉、图像拼接等任务中非常重要。

单应性变换是一种线性变换，用于在不同视角或不同平面之间建立点的对应关系。它可以用一个3×3 的矩阵来表示，称为单应矩阵。这个矩阵描述了从一个图像平面到另一个图像平面的透视变换。

单应性矩阵 (H) 是一个 3×3 的矩阵，用来将一个图像中的点坐标(x,y)映射到另一个图像中的点坐标(x′,y′)。

对点进行归一化处理

import numpy as np

def normalize(points):  #输入的点数组，通常是一个二维数组，其中每一行代表一个点的齐次坐标 ([x, y, w])points=points.astype(float)for row in points:row/=points[-1] #对每个点进行归一化操作。这里的 points[-1] 表示该点的最后一个坐标分量（齐次坐标的权重）。对于每个点，将其所有坐标分量除以该点的权重，以将点归一化到标准形式。return points

points = np.array([[2, 3, 2],[4, 6, 2],[8, 12, 4]])
normalized_points = normalize(points)
print(normalized_points)

在这里插入图片描述

转换为齐次坐标

引入了一个额外的维度来处理点的变换和投影。

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

ones((1,points.shape[1])):创建一个形状为 (1, N) 的数组，其中每个元素都是 1。这一行用于将所有点的齐次坐标的权重分量设置为 1。
np.vstack((points,ones((1,points.shape[1])))):将原始点数组和一行全 1 的数组垂直堆叠，形成一个新的数组，其中每列代表一个点的齐次坐标 [x, y, 1]。

points = np.array([[1, 2, 3],[4, 5, 6]])
homog_points = make_homog(points)
print(homog_points)

在这里插入图片描述

库函数求解单应性矩阵 H

import cv2
import numpy as np
# src_pts 是源图像中的四个点的坐标。
src_pts = np.array([[100, 150], [200, 150], [100, 250], [200, 250]], dtype='float32')
# dst_pts 是目标图像中与 src_pts 对应的四个点的坐标。
dst_pts = np.array([[120, 170], [220, 170], [120, 270], [220, 270]], dtype='float32')
# 这些点用于计算从源图像到目标图像的变换关系。# 计算一个单应性矩阵 H，该矩阵描述了从源图像坐标系到目标图像坐标系的变换。
# status 是一个数组，指示每个点对是否被成功匹配
H, status = cv2.findHomography(src_pts, dst_pts) print("Homography Matrix:\n", H)
print(status)

在这里插入图片描述

验证如下：

在这里插入图片描述

3.1.1 直接线性变换算法（DLT）

用于计算单应性矩阵，基于给定的源点 (fp) 和目标点 (tp)。这个函数主要包含点的归一化、构建方程、求解矩阵以及反归一化的步骤。

from numpy import *

def H_from_points(fp,tp):if fp.shape!=tp.shape: # 确保源点 (fp) 和目标点 (tp) 的形状相同。如果形状不匹配，抛出异常。raise RuntimeError('number of points do not match')# 归一化源点:计算源点的均值m和标准差maxstd。创建归一化矩阵C1，用于将源点fp进行归一化处理，以减小计算中的数值误差。m=mean(fp[:2],axis=1)  # 计算源点的均值 m，对每个坐标分量进行均值计算maxstd=max(std(fp[:2],axis=1))+1e-9    # 计算源点的标准差 maxstd，加一个小偏移量以避免除零错误  C1=diag([1/maxstd,1/maxstd,1])   # 创建归一化矩阵 C1，用于缩放坐标C1[0][2]=-m[0]/maxstd  # 设置 C1 矩阵的平移部分C1[1][2]=-m[1]/maxstd  # 设置 C1 矩阵的平移部分fp=dot(C1,fp)       # 应用归一化矩阵 C1 到源点 fp# 归一化目标点:类似地，对目标点进行归一化处理，计算均值和标准差，并应用归一化矩阵 C2。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)) #构建矩阵A，用于求解单应性矩阵。每对点提供两个方程，总共 2 * nbr_correspondences 行。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) # 使用奇异值分解 (SVD)求解矩阵A的最小特征值对应的向量。取V的最后一行（对应于最小特征值），重塑为3x3矩阵 HH=V[8].reshape((3,3))
#     反归一化H=dot(linalg.inv(C2),dot(H,C1))# 使用逆归一化矩阵将计算得到的单应性矩阵从归一化坐标系转换回原始坐标系，并进行归一化处理。
#     归一化，然后返回return H/H[2,2]

# 定义源点 fp 和目标点 tp
fp = array([[100,200,100,200],[150,150,250,250],[1,1,1,1]])
tp = array([[120,220,120,220],[170,170,270,270],[1,1,1,1]])# 计算单应性矩阵
H = H_from_points(fp, tp)
print(H)
print(dot(H,fp))  #结果为tp

在这里插入图片描述

3.1.2 仿射变换

它是将源点 (fp) 转换到目标点 (tp) 的一种线性变换矩阵。与单应性矩阵不同，仿射变换矩阵不包含透视变换的成分，因此它只适用于进行平移、旋转、缩放和剪切等变换。

def Haffine_from_points(fp,tp):if fp.shape!=tp.shape: # 确保源点 (fp) 和目标点 (tp) 的形状相同。如果形状不匹配，抛出异常。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 矩阵的平移部分C1[1][2]=-m[1]/maxstd  # 设置 C1 矩阵的平移部分fp_cond=dot(C1,fp)  # 应用归一化矩阵 C1 到源点 fpm=mean(tp[:2],axis=1)C2=C1.copy()  #两个点集，必须都进行相同的缩放C2[0][2]=-m[0]/maxstdC2[1][2]=-m[1]/maxstdtp_cond=dot(C2,tp)A=concatenate((fp_cond[:2],tp_cond[:2]),axis=0)  # 拼接源点和目标点的归一化坐标，形成矩阵 AU,S,V=linalg.svd(A.T)tmp=V[:2].T  # 取 V 的前两行，并转置B=tmp[:2]    # 取前两行作为矩阵 BC=tmp[2:4]   # 取接下来的两行作为矩阵 Ctmp2=concatenate((dot(C,linalg.pinv(B)),zeros((2,1))),axis=1) # 计算仿射矩阵的前两列，并添加一列零H=vstack((tmp2,[0,0,1]))   # 形成完整的 3x3 仿射变换矩阵 HH=dot(linalg.inv(C2),dot(H,C1))  # 将仿射矩阵从归一化坐标系转换回原始坐标系return H/H[2,2]

# 定义源点 fp 和目标点 tp
fp = array([[100,200,100,200],[150,150,250,250],[1,1,1,1]])
tp = array([[120,220,120,220],[170,170,270,270],[1,1,1,1]])# 计算仿射变换矩阵
H = Haffine_from_points(fp, tp)
print(H)
print(dot(H,fp))

在这里插入图片描述

3.2 图像扭曲

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

im=array(Image.open('sun.jpg').convert('L'))
H=array([[1.4,0.05,-100],[0.05,1.5,-100],[0,0,1]])#   ndimage.affine_transform()用于应用仿射变换到图像数据上
#   目标将图像im进行坐标变换
#   H[:2,:2]是仿射变换矩阵的前2×2部分,表示旋转和缩放部分，
#   (H[0, 2], H[1, 2])是平移部分，
im2=ndimage.affine_transform(im,H[:2,:2],(H[0,2],H[1,2]))figure(figsize=(10, 3))
gray()
subplot(121)
imshow(im)
subplot(122)
imshow(im2)
show()

在这里插入图片描述

3.2.1 图像中的图像

from numpy import *
from pylab import *
from scipy import ndimage

def image_in_image(im1,im2,tp):m,n=im1.shape[:2] # 获取图像 im1 的高度 (m) 和宽度 (n)fp=array([[0,m,m,0],[0,0,n,n],[1,1,1,1]]) # 定义源图像 im1 的四个角点的齐次坐标
# 创建一个 3x4 的矩阵 fp，它表示源图像 im1 四个角点的齐次坐标。
# 这些点包括左上角 (0,0)，右上角 (m,0)，右下角 (m,n) 和左下角 (0,n)。H=Haffine_from_points(tp,fp)  # 计算从源点 fp 到目标点 tp 的仿射变换矩阵 H。#   ndimage.affine_transform()在原来基础上加上了一个参数
#   im2.shape[:2]:指定了变换后图像的大小。    im1_t=ndimage.affine_transform(im1,H[:2,:2],(H[0,2],H[1,2]),im2.shape[:2])alpha=(im1_t>0)  # 生成一个布尔掩码 alpha，它指示图像 im1_t 中哪些像素值大于 0。这有助于确定哪些区域的图像内容有效。return (1-alpha)*im2+alpha*im1_t
#   将图像 im1_t 和 im2 根据掩码 alpha 进行融合。
#   (1 - alpha) 用于选择目标图像 im2 中的像素，alpha 用于选择变换后的图像 im1_t 中的像素。
#   这样可以将 im1_t 中的有效像素覆盖到 im2 上，并保持 im2 中原有的像素不变。

im1=array(Image.open('beatles.jpg').convert('L'))
im2=array(Image.open('billboard_for_rent.jpg').convert('L'))
figure()
subplot(121)
imshow(im1)
subplot(122)
imshow(im2)# 定义目标点 tp 的坐标，并调用 image_in_image 函数将图像 im2 映射到图像 im1 中，生成变换后的图像 im3
tp=array([[120,260,260,120],   #目标图像中四个点的 y 坐标。[16,16,305,305],    #目标图像中四个点的 x 坐标。[1,1,1,1]])
im3=image_in_image(im1,im2,tp)
# tp 矩阵定义了目标图像的四个角点的位置。仿射变换的目的是将源图像 im1 中的四个角点变换到目标图像中的这些点位置。通过计算这些点之间的仿射变换矩阵，可以将 im1 映射到 im2 中，确保 im1 的四个角点在 im2 中对应到 tp 中指定的位置。figure()
gray()
imshow(im3)
# plot([16,16,305,305],[120,260,260,120],'*')
# axis('equal')
axis('off')
show()

在这里插入图片描述

3.2.2 分段仿射扭曲

from scipy.spatial import Delaunay
from numpy import *
from PIL import Image

def triangulate_points(x, y):""" 二维点的Delaunay三角剖分 """ tri = Delaunay(np.c_[x, y]).simplicesreturn tri

x,y=array(random.standard_normal((2,100)))
# centers,edges,tri,neighbors=Delaunay(x,y)tri = triangulate_points(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()

在这里插入图片描述

def pw_affine(fromim,toim,fp,tp,tri):"""从一幅图像中扭曲矩形图像块fromim = 将要扭曲的图像toim =目标图像fp = 齐次坐标下，扭曲前的点tp =齐次坐标下，扭曲后的点tri = 三角剖分"""im = toim.copy()# 检查图像是灰色图像还是彩色图像is_color = len(fromim.shape) == 3 #创建扭曲后的图像（如果需要对彩色图像的每个颜色通道进行迭代操作，那么有必要这么做）im_t = zeros(im.shape,'uint8')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])#三角形的alphaalpha = alpha_for_triangle(tp[:,t],im.shape[0],im.shape[1])#将三角形加入到图像中im[alpha>0] = im_t[alpha>0]return im

def plot_mesh(x,y,tri):for t in tri:t_ext=[t[0],t[1],t[2],t[0]]plot(x[t_ext],y[t_ext],'r')

fromim=array(Image.open('sunset_tree.jpg'))
x,y=meshgrid(range(5),range(6))
x=(fromim.shape[1]/4)*x.flatten()
y=(fromim.shape[0]/5)*y.flatten()tri=triangulate_points(x,y)im=array(Image.open('turningtorso1.jpg'))
tp=loadtxt('turningtorso1_points.txt')
figure()
subplot(121)
imshow(im)
axis('off')fp=vstack((y,x,ones((1,len(x)))))
tp=vstack((tp[:,1],tp[:,0],ones((1,len(tp)))))im=pw_affine(fromim,im,fp,tp,tri)subplot(122)
imshow(fromim)
for t in tri:t_ext=[t[0],t[1],t[2],t[0]]plot(x[t_ext],y[t_ext],'r')
axis('off')figure()
subplot(121)
imshow(im)
axis('off')subplot(122)
imshow(im)
plot_mesh(tp[1],tp[0],tri)
axis('off')
show()

在这里插入图片描述

3.2.3 图像配准

from xml.dom import minidom
from scipy import linalg
from scipy import ndimage
import os
import imageio
imsave = imageio.imsave

def read_points_from_xml(xmlFileName):""" 读取用于人脸对齐的控制点 """xmldoc = minidom.parse(xmlFileName)facelist = xmldoc.getElementsByTagName('face')faces = {}for xmlFace in facelist:fileName = xmlFace.attributes['file'].valuexf = int(xmlFace.attributes['xf'].value)yf = int(xmlFace.attributes['yf'].value)xs = int(xmlFace.attributes['xs'].value)ys = int(xmlFace.attributes['ys'].value)xm = int(xmlFace.attributes['xm'].value)ym = int(xmlFace.attributes['ym'].value)faces[fileName] = array([xf, yf, xs, ys, xm, ym])return faces

def compute_rigid_transform(refpoints,points):""" 计算用于将点对齐到参考点的旋转、尺度和平移量 """A = array([ [points[0], -points[1], 1, 0],[points[1],  points[0], 0, 1],[points[2], -points[3], 1, 0],[points[3],  points[2], 0, 1],[points[4], -points[5], 1, 0],[points[5],  points[4], 0, 1]])y = array([ refpoints[0],refpoints[1],refpoints[2],refpoints[3],refpoints[4],refpoints[5]])# 计算最小化 ||Ax-y|| 的最小二乘解a,b,tx,ty = linalg.lstsq(A,y)[0]R = array([[a, -b], [b, a]]) # 包含尺度的旋转矩阵return R,tx,ty

def rigid_alignment(faces,path,plotflag=False):""" 严格对齐图像，并将其保存为新的图像path 决定对齐后图像保存的位置设置 plotflag=True，以绘制图像"""# 将第一幅图像中的点作为参考点 
#     print(faces.values())refpoints = list(faces.values())[0]# 使用仿射变换扭曲每幅图像for face in faces:
#         print(os.path.join(path,face))points = faces[face]R,tx,ty = compute_rigid_transform(refpoints, points)T = array([[R[1][1], R[1][0]], [R[0][1], R[0][0]]])im = array(Image.open(os.path.join(path,face)))im2 = zeros(im.shape,'uint8')# 对每个颜色通道进行扭曲for i in range(len(im.shape)):im2[:,:,i] = ndimage.affine_transform(im[:,:,i],linalg.inv(T),offset=[-ty,-tx])if plotflag:imshow(im2)show()# 裁剪边界，并保存对齐后的图像h,w = im2.shape[:2]border = (w+h)/20# 裁剪边界imsave(os.path.join(path,'aligned/'+face),im2[int(border):int(h-border),int(border):int(w-border),:])

# 载入控制点的位置
xmlFileName = r'D:\pyFile\Python计算机视觉编程\data\jkfaces.xml'
points = read_points_from_xml(xmlFileName)
# print(points)
# 注册
rigid_alignment(points,r'D:\pyFile\Python计算机视觉编程\data\jkfaces\\')

figure()
for i in range(1,7):subplot(1,6,i)im=array(Image.open(r'D:\pyFile\Python计算机视觉编程\data\jkfaces\2008010%d.jpg'%i))imshow(im)axis('off')
figure()
for i in range(1,7):subplot(1,6,i)im=array(Image.open(r'D:\pyFile\Python计算机视觉编程\data\jkfaces\aligned\2008010%d.jpg'%i))imshow(im)axis('off')

在这里插入图片描述

3.3 创建全景图

3.3.1 RANSAC

RANSAC是一种用于估计数学模型参数的鲁棒算法，特别是在数据中存在大量异常值时。它最常用于计算机视觉和图像处理中的模型拟合任务。

RANSAC 的工作原理
1. 随机抽样：从数据集中随机选择一小部分数据点，用于计算模型的初始估计。这个小子集通常是模型参数的最小样本量。
2. 模型估计：使用这些随机选择的数据点来拟合模型，计算模型参数。
3. 验证模型：通过将拟合得到的模型应用于所有数据点，确定哪些数据点与模型一致（内点）或不一致（外点）。
4. 评估模型：计算模型的内点数目或模型的其他度量指标。记录模型的内点数目，找出最好的模型。
5. 迭代：重复上述过程多次，每次使用不同的随机样本。选择内点数最多的模型作为最终结果。
6. 优化：一旦确定了最佳模型，使用所有内点来重新估计模型参数，得到更精确的模型。

import numpy
import scipy
import scipy.linalg
def ransac(data,model,n,k,t,d,debug=False,return_all=False):iterations = 0bestfit = Nonebesterr = numpy.infbest_inlier_idxs = Nonewhile iterations < k:maybe_idxs, test_idxs = random_partition(n,data.shape[0])maybeinliers = data[maybe_idxs,:]test_points = data[test_idxs]maybemodel = model.fit(maybeinliers)test_err = model.get_error( test_points, maybemodel)also_idxs = test_idxs[test_err < t] # select indices of rows with accepted pointsalsoinliers = data[also_idxs,:]if debug:print('test_err.min()',test_err.min())print('test_err.max()',test_err.max())print('numpy.mean(test_err)',numpy.mean(test_err))print('iteration %d:len(alsoinliers) = %d'%(iterations,len(alsoinliers)))if len(alsoinliers) > d:betterdata = numpy.concatenate( (maybeinliers, alsoinliers) )bettermodel = model.fit(betterdata)better_errs = model.get_error( betterdata, bettermodel)thiserr = numpy.mean( better_errs )if thiserr < besterr:bestfit = bettermodelbesterr = thiserrbest_inlier_idxs = numpy.concatenate( (maybe_idxs, also_idxs) )iterations+=1if bestfit is None:raise ValueError("did not meet fit acceptance criteria")if return_all:return bestfit, {'inliers':best_inlier_idxs}else:return bestfitdef random_partition(n,n_data):"""return n random rows of data (and also the other len(data)-n rows)"""all_idxs = numpy.arange( n_data )numpy.random.shuffle(all_idxs)idxs1 = all_idxs[:n]idxs2 = all_idxs[n:]return idxs1, idxs2class LinearLeastSquaresModel:"""linear system solved using linear least squaresThis class serves as an example that fulfills the model interfaceneeded by the ransac() function."""def __init__(self,input_columns,output_columns,debug=False):self.input_columns = input_columnsself.output_columns = output_columnsself.debug = debugdef fit(self, data):A = numpy.vstack([data[:,i] for i in self.input_columns]).TB = numpy.vstack([data[:,i] for i in self.output_columns]).Tx,resids,rank,s = numpy.linalg.lstsq(A,B)return xdef get_error( self, data, model):A = numpy.vstack([data[:,i] for i in self.input_columns]).TB = numpy.vstack([data[:,i] for i in self.output_columns]).TB_fit = scipy.dot(A,model)err_per_point = numpy.sum((B-B_fit)**2,axis=1) # sum squared error per rowreturn err_per_pointdef test():# generate perfect input datan_samples = 500n_inputs = 1n_outputs = 1A_exact = 20*numpy.random.random((n_samples,n_inputs) )perfect_fit = 60*numpy.random.normal(size=(n_inputs,n_outputs) ) # the modelB_exact = scipy.dot(A_exact,perfect_fit)assert B_exact.shape == (n_samples,n_outputs)# add a little gaussian noise (linear least squares alone should handle this well)A_noisy = A_exact + numpy.random.normal(size=A_exact.shape )B_noisy = B_exact + numpy.random.normal(size=B_exact.shape )if 1:# add some outliersn_outliers = 100all_idxs = numpy.arange( A_noisy.shape[0] )numpy.random.shuffle(all_idxs)outlier_idxs = all_idxs[:n_outliers]non_outlier_idxs = all_idxs[n_outliers:]A_noisy[outlier_idxs] =  20*numpy.random.random((n_outliers,n_inputs) )B_noisy[outlier_idxs] = 50*numpy.random.normal(size=(n_outliers,n_outputs) )# setup modelall_data = numpy.hstack( (A_noisy,B_noisy) )input_columns = range(n_inputs) # the first columns of the arrayoutput_columns = [n_inputs+i for i in range(n_outputs)] # the last columns of the arraydebug = Truemodel = LinearLeastSquaresModel(input_columns,output_columns,debug=debug)linear_fit,resids,rank,s = numpy.linalg.lstsq(all_data[:,input_columns],all_data[:,output_columns])# run RANSAC algorithmransac_fit, ransac_data = ransac(all_data,model,5, 5000, 7e4, 50, # misc. parametersdebug=debug,return_all=True)if 1:import pylabsort_idxs = numpy.argsort(A_exact[:,0])A_col0_sorted = A_exact[sort_idxs] # maintain as rank-2 arrayif 1:pylab.plot( A_noisy[:,0], B_noisy[:,0], 'k.', label='data' )pylab.plot( A_noisy[ransac_data['inliers'],0], B_noisy[ransac_data['inliers'],0], 'bx', label='RANSAC data' )else:pylab.plot( A_noisy[non_outlier_idxs,0], B_noisy[non_outlier_idxs,0], 'k.', label='noisy data' )pylab.plot( A_noisy[outlier_idxs,0], B_noisy[outlier_idxs,0], 'r.', label='outlier data' )pylab.plot( A_col0_sorted[:,0],numpy.dot(A_col0_sorted,ransac_fit)[:,0],label='RANSAC fit' )pylab.plot( A_col0_sorted[:,0],numpy.dot(A_col0_sorted,perfect_fit)[:,0],label='exact system' )pylab.plot( A_col0_sorted[:,0],numpy.dot(A_col0_sorted,linear_fit)[:,0],label='linear fit' )pylab.legend()pylab.show()if __name__=='__main__':test()

ransac函数：实现了 RANSAC 算法，包含参数设置、模型估计、内点检测等步骤。
LinearLeastSquaresModel类：一个示例模型，使用线性最小二乘法来拟合数据。
test函数：生成测试数据、添加噪声和异常值、运行 RANSAC 算法并可视化结果。

在这里插入图片描述

3.3.2 稳健的单应性矩阵估计

from PIL import Image
from numpy import *
from pylab import *
import os
import subprocess
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['font.sans-serif'] = ['SimHei']  # 黑体字体

以下是第二章SIFT特征匹配所应用的函数，因为本节学习需要，我拷贝了过来，方便学习。（后期是可以封装进行使用的）

def process_image(imagename, resultname, params="--edge-thresh 10 --peak-thresh 5"):if imagename[-3:] != 'pgm':im = Image.open(imagename).convert('L')im.save('tmp.pgm')imagename = 'tmp.pgm'cmmd = str(".\sift.exe " + imagename + " --output=" + resultname + " " + params)os.system(cmmd)print('processed', imagename, 'to', resultname)

def read_features_from_file(filename):f=loadtxt(filename)return f[:,:4],f[:,4:]

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.6desc1_size=desc1.shapematchscores=zeros((desc1_size[0],1),'int')desc2t=desc2.Tfor i in range(desc1_size[0]):dotprods=dot(desc1[i,:],desc2t)dotprods=0.9999*dotprodsindx=argsort(arccos(dotprods))if arccos(dotprods)[indx[0]]<dist_ratio*arccos(dotprods)[indx[1]]:matchscores[i]=int(indx[0])return matchscores

def appendimages(im1,im2):rows1=im1.shape[0]rows2=im2.shape[0]if rows1<rows2:im1=concatenate((im1,zeros((rows2-rows1,im1.shape[1]))),axis=0)elif rows1>rows2:im2=concatenate((im2,zeros((rows1-rows2,im2.shape[1]))),axis=0)return concatenate((im1,im2),axis=1)

def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):print(locs1.shape,locs2.shape)im3=appendimages(im1,im2)  # 将两张图像水平拼接成一张新图像if show_below:im3=vstack((im3,im3))  # 如果 show_below 为 True，将拼接后的图像在垂直方向上再拼接一次figure(figsize=(20, 10))imshow(im3)cols1=im1.shape[1]  # 存储im1的宽度，用于计算绘制线条时的水平偏移量。
#     print(matchscores)for i,m in enumerate(matchscores):  # 会返回一个由索引和值组成的元组value=m[0]if value>0:plot([locs1[i][0],locs2[value][0]+cols1],[locs1[i][1],locs2[value][1]],'c')axis('off')

将这次需要用于全景的五张图片进行处理，提取每张图像的特征，并且匹配相邻图像之间的特征。

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):process_image(imname[i],featname[i])l[i],d[i] = read_features_from_file(featname[i])matches = {}
for i in range(4):matches[i] = match(d[i+1],d[i])

查看图片1和图片2以及匹配结果的可视化。

gray()
im1 = array(Image.open(imname[0]).convert('L'))
im2 = array(Image.open(imname[1]).convert('L'))
plot_matches(im2,im1,l[1],l[0],matches[0])

在这里插入图片描述

查看图片2和图片3以及匹配结果的可视化。

gray()
im1 = array(Image.open(imname[1]).convert('L'))
im2 = array(Image.open(imname[2]).convert('L'))
plot_matches(im2,im1,l[2],l[1],matches[1])

在这里插入图片描述

查看图片3和图片4以及匹配结果的可视化。

gray()
im1 = array(Image.open(imname[3]).convert('L'))
im2 = array(Image.open(imname[2]).convert('L'))
plot_matches(im1,im2,l[3],l[2],matches[2])

在这里插入图片描述

查看图片4和图片5以及匹配结果的可视化。

gray()
im1 = array(Image.open(imname[4]).convert('L'))
im2 = array(Image.open(imname[3]).convert('L'))
plot_matches(im1,im2,l[4],l[3],matches[3])

在这里插入图片描述

接下来使用RANSAC算法求解单应性矩阵H

import numpy as np
import scipy
import scipy.linalg

def H_from_points(fp,tp):if fp.shape!=tp.shape: # 确保源点 (fp) 和目标点 (tp) 的形状相同。如果形状不匹配，抛出异常。raise RuntimeError('number of points do not match')# 归一化源点:计算源点的均值m和标准差maxstd。创建归一化矩阵C1，用于将源点fp进行归一化处理，以减小计算中的数值误差。m=np.mean(fp[:2],axis=1)  # 计算源点的均值 m，对每个坐标分量进行均值计算maxstd=max(np.std(fp[:2],axis=1))+1e-9    # 计算源点的标准差 maxstd，加一个小偏移量以避免除零错误  C1=np.diag([1/maxstd,1/maxstd,1])   # 创建归一化矩阵 C1，用于缩放坐标C1[0][2]=-m[0]/maxstd  # 设置 C1 矩阵的平移部分C1[1][2]=-m[1]/maxstd  # 设置 C1 矩阵的平移部分fp=np.dot(C1,fp)       # 应用归一化矩阵 C1 到源点 fp# 归一化目标点:类似地，对目标点进行归一化处理，计算均值和标准差，并应用归一化矩阵 C2。m=np.mean(tp[:2],axis=1)maxstd=max(np.std(tp[:2],axis=1))+1e-9C2=np.diag([1/maxstd,1/maxstd,1])C2[0][2]=-m[0]/maxstdC2[1][2]=-m[1]/maxstdtp=np.dot(C2,tp)nbr_correspondences=fp.shape[1]A=np.zeros((2*nbr_correspondences,9)) #构建矩阵A，用于求解单应性矩阵。每对点提供两个方程，总共 2 * nbr_correspondences 行。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=np.linalg.svd(A) # 使用奇异值分解 (SVD)求解矩阵A的最小特征值对应的向量。取V的最后一行（对应于最小特征值），重塑为3x3矩阵 HH=V[8].reshape((3,3))
#     反归一化H=np.dot(np.linalg.inv(C2),np.dot(H,C1))# 使用逆归一化矩阵将计算得到的单应性矩阵从归一化坐标系转换回原始坐标系，并进行归一化处理。
#     归一化，然后返回return H/H[2,2]

class RansacModel(object):def __init__(self,debug=False):self.debug=debug def fit(self,data):""" 计算选取的 4 个对应的单应性矩阵 """# 将其转置，来调用 H_from_points() 计算单应性矩阵data = data.T# 映射的起始点fp = data[:3,:4]# 映射的目标点tp = data[3:,:4]# 计算单应性矩阵，然后返回return H_from_points(fp,tp)def get_error(self, data, H):""" 对所有的对应计算单应性矩阵，然后对每个变换后的点，返回相应的误差 """data = data.T# 映射的起始点fp = data[:3]# 映射的目标点tp = data[3:]# 变换fpfp_transformed = dot(H,fp)# 归一化齐次坐标for i in range(3):fp_transformed[i]/=fp_transformed[2]# 返回每个点的误差return sqrt( sum((tp-fp_transformed)**2,axis=0) )

def random_partition(n,n_data):"""return n random rows of data (and also the other len(data)-n rows)"""all_idxs = np.arange( n_data )np.random.shuffle(all_idxs)idxs1 = all_idxs[:n]idxs2 = all_idxs[n:]return idxs1, idxs2

def ransac(data, model, n, k, t, d, debug=False, return_all=False):iterations = 0bestfit = Nonebesterr = np.infbest_inlier_idxs = Nonewhile iterations < k:maybe_idxs, test_idxs = random_partition(n, data.shape[0])maybeinliers = data[maybe_idxs, :]test_points = data[test_idxs]maybemodel = model.fit(maybeinliers)test_err = model.get_error(test_points, maybemodel)also_idxs = test_idxs[test_err < t]  # select indices of rows with accepted pointsalsoinliers = data[also_idxs, :]if debug:print('test_err.min()', test_err.min())print('test_err.max()', test_err.max())print('numpy.mean(test_err)', np.mean(test_err))print('iteration %d:len(alsoinliers) = %d' %(iterations, len(alsoinliers)))if len(alsoinliers) > d:betterdata = np.concatenate((maybeinliers, alsoinliers))bettermodel = model.fit(betterdata)better_errs = model.get_error(betterdata, bettermodel)#重新计算总的errorthiserr = np.mean(better_errs)if thiserr < besterr:bestfit = bettermodelbesterr = thiserrbest_inlier_idxs = np.concatenate((maybe_idxs, also_idxs))iterations += 1if bestfit is None:raise ValueError("did not meet fit acceptance criteria")if return_all:return bestfit, {'inliers': best_inlier_idxs}else:return bestfit

def H_from_ransac(fp,tp,model,maxiter=1000,match_theshold=10):# 对应点组data = vstack((fp,tp))# 计算 H，并返回H,ransac_data = ransac(data.T,model,4,maxiter,match_theshold,10,return_all=True)
#     print(H,ransac_data)return H,ransac_data['inliers']

def make_homog(points):  # 输入的二维点数组，通常形状为(2,N)，其中N是点的数量。每一列代表一个二维点的坐标[x,y]return np.vstack((points,np.ones((1,points.shape[1]))))

# 将匹配转换成齐次坐标点的函数
def convert_points(j):ndx = matches[j].nonzero()[0]fp = make_homog(l[j+1][ndx,:2].T)ndx2 = [int(matches[j][i]) for i in ndx]tp = make_homog(l[j][ndx2,:2].T)return fp,tp# 估计单应性矩阵
model = RansacModel()fp,tp = convert_points(1)
H_12 = H_from_ransac(fp,tp,model)[0] # im1 到 im2 的单应性矩阵
# print(tp)
# print(dot(H_12,fp))
# print('---')fp,tp = convert_points(0)
H_01 = H_from_ransac(fp,tp,model)[0] # im0 到 im1 的单应性矩阵tp,fp = convert_points(2) # 注意:点是反序的
H_32 = H_from_ransac(fp,tp,model)[0] # im3 到 im2 的单应性矩阵tp,fp = convert_points(3) # 注意:点是反序的
H_43 = H_from_ransac(fp,tp,model)[0] # im4 到 im3 的单应性矩阵

3.3.3 拼接图像

from scipy import ndimage

接下来是panorama全景图的主要函数，这里需要注意的是，或许是因为我使用的是jupyter，所有x和y的坐标有些许相反，我按照课本上的函数执行是会发生位置对调（也是调了半天），所以我把三处和课本不一致的地方标了出来

def panorama(H,fromim,toim,padding=2400,delta=2400):is_color = len(fromim.shape) == 3  #用于检查图像是灰度图像还是彩色图像（3通道）# 用于geometric_transform()的单应性变换def transf(p):  # 输入点的坐标，通常是一个二维坐标 (x, y)，表示图像上的一个点。p2 = dot(H,[p[1],p[0],1])  # p[1]为x坐标，p[0]为y坐标，跟课本上的第一处不一样) 
#         p2 = dot(H,[p[0],p[1],1])
#         print(p)# 与上面相同，进行归一化后，p2[1]/p2[2]为x坐标，p2[0]/p2[2]为y坐标，跟课本上的第二处不一样) return (p2[1]/p2[2],p2[0]/p2[2])
#         return (p2[1]/p2[2],p2[0]/p2[2])if H[0,2]<0: # fromim在右边,H[0,2]为y轴的偏移量(跟课本上的第三处不一样) 
#   if H[1,2]<0:print('warp - right')# 变换fromimif is_color:# 在目标图像的右边填充0# hstack:水平拼接数组toim_t=hstack((toim,zeros((toim.shape[0],padding,3))))#  fromim_t用于存储变换后的 fromim 图像。这个数组的宽度是目标图像宽度加上 padding，以容纳图像拼接后的结果。fromim_t=zeros((toim.shape[0],toim.shape[1]+padding,toim.shape[2]))for col in range(3):
#               ndimage.geometric_transform 会对图像中的每一个像素应用变换函数 transf。fromim_t[:,:,col]=ndimage.geometric_transform(fromim[:,:,col],transf,(toim.shape[0],toim.shape[1]+padding))else:# 在目标图像的右边填充0toim_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_delta = array([[1,0,-delta],[0,1,0],[0,0,1]])H = dot(H,H_delta)# fromim变换if is_color:# 在目标图像的左边填充0toim_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:# 在目标图像的左边填充0toim_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

以第3张图为中心，图1和图2拼接在右边，图4和图5拼接在左边

# 扭曲图像
delta = 2000 # 用于填充和平移
im1 = array(Image.open(imname[1]))
im2 = array(Image.open(imname[2]))
im_12 = panorama(H_12,im1,im2,delta,delta)im1 = array(Image.open(imname[0]))
im_02 = panorama(dot(H_12,H_01),im1,im_12,delta,delta)im1 = array(Image.open(imname[3]))
im_32 = panorama(H_32,im1,im_02,delta,delta)im1 = array(Image.open(imname[4]))
im_42 = panorama(dot(H_32,H_43),im1,im_32,delta,2*delta)imshow(array(im_42, "uint8"))
axis('off')
show()