滨州做网站建设的公司,网络营销上市公司,视频网站如何做推广,wordpress添加自定义tag标签python计算机视觉编程——3.图像到图像的映射 3.图像到图像的映射3.1 单应性变换3.1.1 直接线性变换算法#xff08;DLT#xff09;3.1.2 仿射变换 3.2 图像扭曲3.2.1 图像中的图像3.2.2 分段仿射扭曲3.2.3 图像配准 3.3 创建全景图3.3.1 RANSAC3.3.2 稳健的单应性矩阵估计3.… python计算机视觉编程——3.图像到图像的映射 3.图像到图像的映射3.1 单应性变换3.1.1 直接线性变换算法DLT3.1.2 仿射变换 3.2 图像扭曲3.2.1 图像中的图像3.2.2 分段仿射扭曲3.2.3 图像配准 3.3 创建全景图3.3.1 RANSAC3.3.2 稳健的单应性矩阵估计3.3.3 拼接图像 3.图像到图像的映射
3.1 单应性变换
将一个图像平面上的点映射到另一个图像平面上的技术。它在图像配准、立体视觉、图像拼接等任务中非常重要。
单应性变换是一种线性变换用于在不同视角或不同平面之间建立点的对应关系。它可以用一个3×3 的矩阵来表示称为单应矩阵。这个矩阵描述了从一个图像平面到另一个图像平面的透视变换。
单应性矩阵 (H) 是一个 3×3 的矩阵用来将一个图像中的点坐标(x,y)映射到另一个图像中的点坐标(x′,y′)。
对点进行归一化处理
import numpy as npdef normalize(points): #输入的点数组通常是一个二维数组其中每一行代表一个点的齐次坐标 ([x, y, w])pointspoints.astype(float)for row in points:row/points[-1] #对每个点进行归一化操作。这里的 points[-1] 表示该点的最后一个坐标分量齐次坐标的权重。对于每个点将其所有坐标分量除以该点的权重以将点归一化到标准形式。return pointspoints 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]], dtypefloat32)
# dst_pts 是目标图像中与 src_pts 对应的四个点的坐标。
dst_pts np.array([[120, 170], [220, 170], [120, 270], [220, 270]], dtypefloat32)
# 这些点用于计算从源图像到目标图像的变换关系。# 计算一个单应性矩阵 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进行归一化处理以减小计算中的数值误差。mmean(fp[:2],axis1) # 计算源点的均值 m对每个坐标分量进行均值计算maxstdmax(std(fp[:2],axis1))1e-9 # 计算源点的标准差 maxstd加一个小偏移量以避免除零错误 C1diag([1/maxstd,1/maxstd,1]) # 创建归一化矩阵 C1用于缩放坐标C1[0][2]-m[0]/maxstd # 设置 C1 矩阵的平移部分C1[1][2]-m[1]/maxstd # 设置 C1 矩阵的平移部分fpdot(C1,fp) # 应用归一化矩阵 C1 到源点 fp# 归一化目标点:类似地对目标点进行归一化处理计算均值和标准差并应用归一化矩阵 C2。mmean(tp[:2],axis1)maxstdmax(std(tp[:2],axis1))1e-9C2diag([1/maxstd,1/maxstd,1])C2[0][2]-m[0]/maxstdC2[1][2]-m[1]/maxstdtpdot(C2,tp)nbr_correspondencesfp.shape[1]Azeros((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*i1][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,Vlinalg.svd(A) # 使用奇异值分解 (SVD)求解矩阵A的最小特征值对应的向量。取V的最后一行对应于最小特征值重塑为3x3矩阵 HHV[8].reshape((3,3))
# 反归一化Hdot(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)) #结果为tp3.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)mmean(fp[:2],axis1)maxstdmax(std(fp[:2],axis1))1e-9 C1diag([1/maxstd,1/maxstd,1])C1[0][2]-m[0]/maxstd # 设置 C1 矩阵的平移部分C1[1][2]-m[1]/maxstd # 设置 C1 矩阵的平移部分fp_conddot(C1,fp) # 应用归一化矩阵 C1 到源点 fpmmean(tp[:2],axis1)C2C1.copy() #两个点集必须都进行相同的缩放C2[0][2]-m[0]/maxstdC2[1][2]-m[1]/maxstdtp_conddot(C2,tp)Aconcatenate((fp_cond[:2],tp_cond[:2]),axis0) # 拼接源点和目标点的归一化坐标形成矩阵 AU,S,Vlinalg.svd(A.T)tmpV[:2].T # 取 V 的前两行并转置Btmp[:2] # 取前两行作为矩阵 BCtmp[2:4] # 取接下来的两行作为矩阵 Ctmp2concatenate((dot(C,linalg.pinv(B)),zeros((2,1))),axis1) # 计算仿射矩阵的前两列并添加一列零Hvstack((tmp2,[0,0,1])) # 形成完整的 3x3 仿射变换矩阵 HHdot(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 *imarray(Image.open(sun.jpg).convert(L))
Harray([[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])是平移部分
im2ndimage.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 ndimagedef image_in_image(im1,im2,tp):m,nim1.shape[:2] # 获取图像 im1 的高度 (m) 和宽度 (n)fparray([[0,m,m,0],[0,0,n,n],[1,1,1,1]]) # 定义源图像 im1 的四个角点的齐次坐标
# 创建一个 3x4 的矩阵 fp它表示源图像 im1 四个角点的齐次坐标。
# 这些点包括左上角 (0,0)右上角 (m,0)右下角 (m,n) 和左下角 (0,n)。HHaffine_from_points(tp,fp) # 计算从源点 fp 到目标点 tp 的仿射变换矩阵 H。# ndimage.affine_transform()在原来基础上加上了一个参数
# im2.shape[:2]:指定了变换后图像的大小。 im1_tndimage.affine_transform(im1,H[:2,:2],(H[0,2],H[1,2]),im2.shape[:2])alpha(im1_t0) # 生成一个布尔掩码 alpha它指示图像 im1_t 中哪些像素值大于 0。这有助于确定哪些区域的图像内容有效。return (1-alpha)*im2alpha*im1_t
# 将图像 im1_t 和 im2 根据掩码 alpha 进行融合。
# (1 - alpha) 用于选择目标图像 im2 中的像素alpha 用于选择变换后的图像 im1_t 中的像素。
# 这样可以将 im1_t 中的有效像素覆盖到 im2 上并保持 im2 中原有的像素不变。im1array(Image.open(beatles.jpg).convert(L))
im2array(Image.open(billboard_for_rent.jpg).convert(L))
figure()
subplot(121)
imshow(im1)
subplot(122)
imshow(im2)# 定义目标点 tp 的坐标并调用 image_in_image 函数将图像 im2 映射到图像 im1 中生成变换后的图像 im3
tparray([[120,260,260,120], #目标图像中四个点的 y 坐标。[16,16,305,305], #目标图像中四个点的 x 坐标。[1,1,1,1]])
im3image_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 Imagedef triangulate_points(x, y): 二维点的Delaunay三角剖分 tri Delaunay(np.c_[x, y]).simplicesreturn trix,yarray(random.standard_normal((2,100)))
# centers,edges,tri,neighborsDelaunay(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[alpha0] im_t[alpha0]return imdef 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)fromimarray(Image.open(sunset_tree.jpg))
x,ymeshgrid(range(5),range(6))
x(fromim.shape[1]/4)*x.flatten()
y(fromim.shape[0]/5)*y.flatten()tritriangulate_points(x,y)imarray(Image.open(turningtorso1.jpg))
tploadtxt(turningtorso1_points.txt)
figure()
subplot(121)
imshow(im)
axis(off)fpvstack((y,x,ones((1,len(x)))))
tpvstack((tp[:,1],tp[:,0],ones((1,len(tp)))))impw_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.imsavedef 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 facesdef 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,tydef rigid_alignment(faces,path,plotflagFalse): 严格对齐图像并将其保存为新的图像path 决定对齐后图像保存的位置设置 plotflagTrue以绘制图像# 将第一幅图像中的点作为参考点
# 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 (wh)/20# 裁剪边界imsave(os.path.join(path,aligned/face),im2[int(border):int(h-border),int(border):int(w-border),:])# 载入控制点的位置
xmlFileName rD:\pyFile\Python计算机视觉编程\data\jkfaces.xml
points read_points_from_xml(xmlFileName)
# print(points)
# 注册
rigid_alignment(points,rD:\pyFile\Python计算机视觉编程\data\jkfaces\\)figure()
for i in range(1,7):subplot(1,6,i)imarray(Image.open(rD:\pyFile\Python计算机视觉编程\data\jkfaces\2008010%d.jpg%i))imshow(im)axis(off)
figure()
for i in range(1,7):subplot(1,6,i)imarray(Image.open(rD:\pyFile\Python计算机视觉编程\data\jkfaces\aligned\2008010%d.jpg%i))imshow(im)axis(off)3.3 创建全景图
3.3.1 RANSAC
RANSAC是一种用于估计数学模型参数的鲁棒算法特别是在数据中存在大量异常值时。它最常用于计算机视觉和图像处理中的模型拟合任务。
RANSAC 的工作原理 随机抽样 从数据集中随机选择一小部分数据点用于计算模型的初始估计。这个小子集通常是模型参数的最小样本量。模型估计 使用这些随机选择的数据点来拟合模型计算模型参数。验证模型 通过将拟合得到的模型应用于所有数据点确定哪些数据点与模型一致内点或不一致外点。评估模型 计算模型的内点数目或模型的其他度量指标。记录模型的内点数目找出最好的模型。迭代 重复上述过程多次每次使用不同的随机样本。选择内点数最多的模型作为最终结果。优化 一旦确定了最佳模型使用所有内点来重新估计模型参数得到更精确的模型。
import numpy
import scipy
import scipy.linalg
def ransac(data,model,n,k,t,d,debugFalse,return_allFalse):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) )iterations1if 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,debugFalse):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,axis1) # 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(sizeA_exact.shape )B_noisy B_exact numpy.random.normal(sizeB_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_inputsi for i in range(n_outputs)] # the last columns of the arraydebug Truemodel LinearLeastSquaresModel(input_columns,output_columns,debugdebug)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. parametersdebugdebug,return_allTrue)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., labeldata )pylab.plot( A_noisy[ransac_data[inliers],0], B_noisy[ransac_data[inliers],0], bx, labelRANSAC data )else:pylab.plot( A_noisy[non_outlier_idxs,0], B_noisy[non_outlier_idxs,0], k., labelnoisy data )pylab.plot( A_noisy[outlier_idxs,0], B_noisy[outlier_idxs,0], r., labeloutlier data )pylab.plot( A_col0_sorted[:,0],numpy.dot(A_col0_sorted,ransac_fit)[:,0],labelRANSAC fit )pylab.plot( A_col0_sorted[:,0],numpy.dot(A_col0_sorted,perfect_fit)[:,0],labelexact system )pylab.plot( A_col0_sorted[:,0],numpy.dot(A_col0_sorted,linear_fit)[:,0],labellinear 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.pgmcmmd str(.\sift.exe imagename --output resultname params)os.system(cmmd)print(processed, imagename, to, resultname)def read_features_from_file(filename):floadtxt(filename)return f[:,:4],f[:,4:]def match(desc1,desc2):desc1array([d/linalg.norm(d) for d in desc1])desc2array([d/linalg.norm(d) for d in desc2])dist_ratio0.6desc1_sizedesc1.shapematchscoreszeros((desc1_size[0],1),int)desc2tdesc2.Tfor i in range(desc1_size[0]):dotprodsdot(desc1[i,:],desc2t)dotprods0.9999*dotprodsindxargsort(arccos(dotprods))if arccos(dotprods)[indx[0]]dist_ratio*arccos(dotprods)[indx[1]]:matchscores[i]int(indx[0])return matchscoresdef appendimages(im1,im2):rows1im1.shape[0]rows2im2.shape[0]if rows1rows2:im1concatenate((im1,zeros((rows2-rows1,im1.shape[1]))),axis0)elif rows1rows2:im2concatenate((im2,zeros((rows1-rows2,im2.shape[1]))),axis0)return concatenate((im1,im2),axis1)def plot_matches(im1,im2,locs1,locs2,matchscores,show_belowTrue):print(locs1.shape,locs2.shape)im3appendimages(im1,im2) # 将两张图像水平拼接成一张新图像if show_below:im3vstack((im3,im3)) # 如果 show_below 为 True将拼接后的图像在垂直方向上再拼接一次figure(figsize(20, 10))imshow(im3)cols1im1.shape[1] # 存储im1的宽度用于计算绘制线条时的水平偏移量。
# print(matchscores)for i,m in enumerate(matchscores): # 会返回一个由索引和值组成的元组valuem[0]if value0:plot([locs1[i][0],locs2[value][0]cols1],[locs1[i][1],locs2[value][1]],c)axis(off)将这次需要用于全景的五张图片进行处理提取每张图像的特征并且匹配相邻图像之间的特征。
featname [Univstr(i1).sift for i in range(5)]
imname [Univstr(i1).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[i1],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.linalgdef H_from_points(fp,tp):if fp.shape!tp.shape: # 确保源点 (fp) 和目标点 (tp) 的形状相同。如果形状不匹配抛出异常。raise RuntimeError(number of points do not match)# 归一化源点:计算源点的均值m和标准差maxstd。创建归一化矩阵C1用于将源点fp进行归一化处理以减小计算中的数值误差。mnp.mean(fp[:2],axis1) # 计算源点的均值 m对每个坐标分量进行均值计算maxstdmax(np.std(fp[:2],axis1))1e-9 # 计算源点的标准差 maxstd加一个小偏移量以避免除零错误 C1np.diag([1/maxstd,1/maxstd,1]) # 创建归一化矩阵 C1用于缩放坐标C1[0][2]-m[0]/maxstd # 设置 C1 矩阵的平移部分C1[1][2]-m[1]/maxstd # 设置 C1 矩阵的平移部分fpnp.dot(C1,fp) # 应用归一化矩阵 C1 到源点 fp# 归一化目标点:类似地对目标点进行归一化处理计算均值和标准差并应用归一化矩阵 C2。mnp.mean(tp[:2],axis1)maxstdmax(np.std(tp[:2],axis1))1e-9C2np.diag([1/maxstd,1/maxstd,1])C2[0][2]-m[0]/maxstdC2[1][2]-m[1]/maxstdtpnp.dot(C2,tp)nbr_correspondencesfp.shape[1]Anp.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*i1][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,Vnp.linalg.svd(A) # 使用奇异值分解 (SVD)求解矩阵A的最小特征值对应的向量。取V的最后一行对应于最小特征值重塑为3x3矩阵 HHV[8].reshape((3,3))
# 反归一化Hnp.dot(np.linalg.inv(C2),np.dot(H,C1))# 使用逆归一化矩阵将计算得到的单应性矩阵从归一化坐标系转换回原始坐标系并进行归一化处理。
# 归一化然后返回return H/H[2,2]class RansacModel(object):def __init__(self,debugFalse):self.debugdebug 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,axis0) )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, idxs2def ransac(data, model, n, k, t, d, debugFalse, return_allFalse):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 bestfitdef H_from_ransac(fp,tp,model,maxiter1000,match_theshold10):# 对应点组data vstack((fp,tp))# 计算 H并返回H,ransac_data ransac(data.T,model,4,maxiter,match_theshold,10,return_allTrue)
# 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[j1][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,padding2400,delta2400):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_thstack((toim,zeros((toim.shape[0],padding,3))))# fromim_t用于存储变换后的 fromim 图像。这个数组的宽度是目标图像宽度加上 padding以容纳图像拼接后的结果。fromim_tzeros((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()
