python计算机视觉编程——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-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)) #构建矩阵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矩阵 H
H=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 到源点 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)
A=concatenate((fp_cond[:2],tp_cond[:2]),axis=0) # 拼接源点和目标点的归一化坐标,形成矩阵 A
U,S,V=linalg.svd(A.T)
tmp=V[:2].T # 取 V 的前两行,并转置
B=tmp[:2] # 取前两行作为矩阵 B
C=tmp[2:4] # 取接下来的两行作为矩阵 C
tmp2=concatenate((dot(C,linalg.pinv(B)),zeros((2,1))),axis=1) # 计算仿射矩阵的前两列,并添加一列零
H=vstack((tmp2,[0,0,1])) # 形成完整的 3x3 仿射变换矩阵 H
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 = 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]).simplices
return 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])
#三角形的alpha
alpha = 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'].value
xf = 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 的工作原理
- 随机抽样: 从数据集中随机选择一小部分数据点,用于计算模型的初始估计。这个小子集通常是模型参数的最小样本量。
- 模型估计: 使用这些随机选择的数据点来拟合模型,计算模型参数。
- 验证模型: 通过将拟合得到的模型应用于所有数据点,确定哪些数据点与模型一致(内点)或不一致(外点)。
- 评估模型: 计算模型的内点数目或模型的其他度量指标。记录模型的内点数目,找出最好的模型。
- 迭代: 重复上述过程多次,每次使用不同的随机样本。选择内点数最多的模型作为最终结果。
- 优化: 一旦确定了最佳模型,使用所有内点来重新估计模型参数,得到更精确的模型。
import numpy
import scipy
import scipy.linalg
def ransac(data,model,n,k,t,d,debug=False,return_all=False):
iterations = 0
bestfit = None
besterr = numpy.inf
best_inlier_idxs = None
while 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 points
alsoinliers = 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 = bettermodel
besterr = thiserr
best_inlier_idxs = numpy.concatenate( (maybe_idxs, also_idxs) )
iterations+=1
if bestfit is None:
raise ValueError("did not meet fit acceptance criteria")
if return_all:
return bestfit, {'inliers':best_inlier_idxs}
else:
return bestfit
def 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, idxs2
class LinearLeastSquaresModel:
"""linear system solved using linear least squares
This class serves as an example that fulfills the model interface
needed by the ransac() function.
"""
def __init__(self,input_columns,output_columns,debug=False):
self.input_columns = input_columns
self.output_columns = output_columns
self.debug = debug
def fit(self, data):
A = numpy.vstack([data[:,i] for i in self.input_columns]).T
B = numpy.vstack([data[:,i] for i in self.output_columns]).T
x,resids,rank,s = numpy.linalg.lstsq(A,B)
return x
def get_error( self, data, model):
A = numpy.vstack([data[:,i] for i in self.input_columns]).T
B = numpy.vstack([data[:,i] for i in self.output_columns]).T
B_fit = scipy.dot(A,model)
err_per_point = numpy.sum((B-B_fit)**2,axis=1) # sum squared error per row
return err_per_point
def test():
# generate perfect input data
n_samples = 500
n_inputs = 1
n_outputs = 1
A_exact = 20*numpy.random.random((n_samples,n_inputs) )
perfect_fit = 60*numpy.random.normal(size=(n_inputs,n_outputs) ) # the model
B_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 outliers
n_outliers = 100
all_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 model
all_data = numpy.hstack( (A_noisy,B_noisy) )
input_columns = range(n_inputs) # the first columns of the array
output_columns = [n_inputs+i for i in range(n_outputs)] # the last columns of the array
debug = True
model = 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 algorithm
ransac_fit, ransac_data = ransac(all_data,model,
5, 5000, 7e4, 50, # misc. parameters
debug=debug,return_all=True)
if 1:
import pylab
sort_idxs = numpy.argsort(A_exact[:,0])
A_col0_sorted = A_exact[sort_idxs] # maintain as rank-2 array
if 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.6
desc1_size=desc1.shape
matchscores=zeros((desc1_size[0],1),'int')
desc2t=desc2.T
for i in range(desc1_size[0]):
dotprods=dot(desc1[i,:],desc2t)
dotprods=0.9999*dotprods
indx=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-9
C2=np.diag([1/maxstd,1/maxstd,1])
C2[0][2]=-m[0]/maxstd
C2[1][2]=-m[1]/maxstd
tp=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矩阵 H
H=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:]
# 变换fp
fp_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 = 0
bestfit = None
besterr = np.inf
best_inlier_idxs = None
while 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 points
alsoinliers = 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)#重新计算总的error
thiserr = np.mean(better_errs)
if thiserr < besterr:
bestfit = bettermodel
besterr = thiserr
best_inlier_idxs = np.concatenate((maybe_idxs, also_idxs))
iterations += 1
if 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')
# 变换fromim
if 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:
# 在目标图像的右边填充0
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_delta = array([[1,0,-delta],[0,1,0],[0,0,1]])
H = dot(H,H_delta)
# fromim变换
if is_color:
# 在目标图像的左边填充0
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:
# 在目标图像的左边填充0
toim_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()