python计算机视觉编程第三章——图像到图像的映射

小小怪下士(๑◝ᴗ◜๑)

已于 2024-09-03 15:06:12 修改

阅读量399

点赞数 4

文章标签： python 计算机视觉开发语言

于 2024-08-27 13:42:00 首次发布

本文链接：https://blog.csdn.net/qq_67725022/article/details/141563201

版权

一、单应性变换

单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。在这里，平面是指图像或者三维中的平面表面。单应性变换具有很强的实用性，比如图像配准、图像纠正和纹理扭曲，以及创建全景图像。我们将频繁地使用单应性变换。

下面的函数可以实现对点进行归一化和转换齐次坐标的功能:

def normalize(points):
    """ 在齐次坐标意义下，对点集进行归一化，使最后一行为 1 """
    for row in points:
        row /= points[-1]
    return points

def make_homog(points):
    """ 将点集(dim×n 的数组)转换为齐次坐标表示 """

    return vstack((points,ones((1,points.shape[1]))))

1直接线性变换算法

DLT(Direct Linear Transformation，直接线性变换)是给定4个或者更多对应点对矩阵，来计算单应性矩阵 H 的算法。将单应性矩阵 H 作用在对应点对上，重新写出该方程，我们可以得到下面的方程:

或者 Ah=0，其中 A 是一个具有对应点对二倍数量行数的矩阵。将这些对应点对方程的系数堆叠到一个矩阵中，我们可以使用 SVD(Singular Value Decomposition，奇异值分解)算法找到 H 的最小二乘解。下面是该算法的代码:

def H_from_points(fp,tp):
    """ 使用线性 DLT 方法，计算单应性矩阵 H，使 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-9 C1 = diag([1/maxstd, 1/maxstd, 1]) C1[0][2] = -m[0]/maxstd
    C1[1][2] = -m[1]/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[1]/maxstd
    tp = 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[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]

2仿射变换

由于仿射变换具有 6 个自由度，因此我们需要三个对应点对来估计矩阵 H。通过将最后两个元素设置为 0，即 h7=h8=0，仿射变换可以用上面的 DLT 算法估计得出。

def Haffine_from_points(fp,tp):
    """ 计算 H，仿射变换，使得 tp 是 fp 经过仿射变换 H 得到的 """
    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[1]/maxstd
        fp_cond = dot(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)

        # 因为归一化后点的均值为 0，所以平移量为 0
        A = concatenate((fp_cond[:2],tp_cond[:2]), axis=0)
        U,S,V = linalg.svd(A.T)

        # 如 Hartley 和 Zisserman 著的 Multiple View Geometry in Computer , Scond Edition 所示，
        # 创建矩阵B和C
        tmp = V[:2].T
        B = 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]

二、图像扭曲

from scipy import ndimage
from PIL import Image
from numpy import *
from matplotlib.pyplot import *
from matplotlib import pyplot as plt

im = array(Image.open('../img/eye.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]))
figure()
gray()
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.subplot(121), plt.imshow(im, plt.cm.gray), plt.title('原始灰度图像'), plt.axis('off')
plt.subplot(122), plt.imshow(im2, plt.cm.gray), plt.title('扭曲后图像'), plt.axis('off')
show()

1图像中的图像

def image_in_image(im1,im2,tp):
 """ 使用仿射变换将 im1 放置在 im2 上，使 im1 图像的角和 tp 尽可能的靠近
 tp 是齐次表示的，并且是按照从左上角逆时针计算的 """
 # 扭曲的点
 m,n = im1.shape[:2]
 fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])
 # 计算仿射变换，并且将其应用于图像 im1
 H = homography.Haffine_from_points(tp,fp)
 im1_t = ndimage.affine_transform(im1,H[:2,:2],
 (H[0,2],H[1,2]),im2.shape[:2])
 alpha = (im1_t > 0)
 return (1-alpha)*im2 + alpha*im1_t

# 仿射扭曲 im1 到 im2 的例子
im1 = array(Image.open('../img/eye.jpg').convert('L'))
im2 = array(Image.open('../img/crans_1_small.jpg').convert('L'))
# 选定一些目标点
tp3 = array([[364,640,640,364],[40,40,605,605],[1,1,1,1]])
im3 = warp.image_in_image(im1,im2,tp3)
tp4 = array([[64,340,340,64],[340,340,905,905],[1,1,1,1]])
im4 = warp.image_in_image(im1,im2,tp4)
figure()
gray()
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.subplot(121), plt.imshow(im3, plt.cm.gray), plt.axis('off')
plt.subplot(122), plt.imshow(im4, plt.cm.gray), plt.axis('off')
show()

import warp

# 仿射扭曲 im1 到 im2 的例子
from py_cvp import homography

im1 = array(Image.open('../img/eye.jpg').convert('L'))
im2 = array(Image.open('../img/ggp.jpg').convert('L'))

tp = array([[30,140,125,4],[130,120,430,430],[1,1,1,1]])

im3 = warp.image_in_image(im1,im2,tp)

# 选定 im1 角上的一些点
m,n = im1.shape[:2]
fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])
# 第一个三角形
tp2 = tp[:,:3]
fp2 = fp[:,:3]
# 计算 H
H = homography.Haffine_from_points(tp2,fp2)
im1_t = ndimage.affine_transform(im1,H[:2,:2],
 (H[0,2],H[1,2]),im2.shape[:2])
# 三角形的 alpha
alpha = warp.alpha_for_triangle(tp2,im2.shape[0],im2.shape[1])
im4_1 = (1-alpha)*im2 + alpha*im1_t
# 第二个三角形
tp2 = tp[:,[0,2,3]]
fp2 = fp[:,[0,2,3]]
# 计算 H
H = homography.Haffine_from_points(tp2,fp2)
im1_t = ndimage.affine_transform(im1,H[:2,:2],
 (H[0,2],H[1,2]),im2.shape[:2])
# 三角形的 alpha 图像
alpha = warp.alpha_for_triangle(tp2,im2.shape[0],im2.shape[1])
im4 = (1-alpha)*im4_1 + alpha*im1_t

figure()
gray()
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.subplot(121), plt.imshow(im3, plt.cm.gray),plt.title('仿射扭曲'),plt.axis('off')
plt.subplot(122), plt.imshow(im4, plt.cm.gray),plt.title('使用两个三角形的仿射弯曲'), plt.axis('off')
show()

2分段仿射扭曲

三角形图像块的仿射扭曲可以完成角点的精确匹配。对应点对集合之间最常用的扭曲方式是分段仿射扭曲。可以使用使用 Matplotlib 和 SciPy 来完成这个操作。为了三角化这些点，经常使用狄洛克三角剖分方法。其一个简单的例子为：

运行错误：

import matplotlib.delaunay as md

报错：ModuleNotFoundError: No module named ‘matplotlib.delaunay’

改为：from scipy.spatial import Delaunay

并把triangulate_points(x,y)里面的代码替换成
tri = Delaunay(np.c_[x,y]).simplices

使用狄洛克三角剖分标记点进行分段仿射扭曲：

主要代码：

from scipy.spatial import Delaunay

def triangulate_points(x,y):
 """ 二维点的 Delaunay 三角剖分 """
 #centers,edges,tri,neighbors = md.delaunay(x,y)
 tri = Delaunay(np.c_[x, y]).simplices
 return tri



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 = homography.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')

使用ginput()手工取点：

运行错误：

‘numpy.float64‘ object cannot be interpreted as an integer

 for i in range(min(points[0]),max(points[0])):
  for j in range(min(points[1]),max(points[1])):

改为 for i in range(min(points[0]).astype(int),max(points[0]).astype(int)):
  for j in range(min(points[1]).astype(int),max(points[1]).astype(int)):

3图像配准

主要代码：

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




from scipy import linalg
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，以绘制图像 """
 # 将第一幅图像中的点作为参考点
 refpoints = list(faces.values())[0]
 print(refpoints)
 # 使用仿射变换扭曲每幅图像
 for face in faces:
    points = faces[face]
    R,tx,ty = compute_rigid_transform(refpoints, points)
    T = ar

三、创建全景图

1RANSAC

RANSAC是“RANdom SAmple Consensus”(随机一致性采样）的缩写。该方法是用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型，例如点集之间的单应性矩阵，RANSAC基本的思想是，数据中包含正确的点和噪声点，合理的模型应该能够在描述正确数据点的同时摒弃噪声点。

2稳健的单应性矩阵估计

为使用 SIFT 特征自动找到匹配对应：

import sift
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):
    sift.process_image(imname[i],featname[i])
    l[i],d[i] = sift.read_features_from_file(featname[i])

matches = {}
for i in range(4):
    matches[i] = sift.match(d[i+1],d[i])

关键代码：

class RansacModel(object):
    """ 用于测试单应性矩阵的类，其中单应性矩阵是由网站 http://www.scipy.org/Cookbook/RANSAC 上
    的 ransac.py 计算出来的
    """

    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 H_from_ransac(fp,tp,model,maxiter=1000,match_theshold=10):
    """ 使用 RANSAC 稳健性估计点对应间的单应性矩阵 H(ransac.py 为从
        http://www.scipy.org/Cookbook/RANSAC 下载的版本)
        # 输入:齐次坐标表示的点 fp，tp(3×n 的数组)"""

        import ransac

        # 对应点组
        data = vstack((fp,tp))
        # 计算 H，并返回
        H,ransac_data = ransac.ransac(data.T,model,4,maxiter,match_theshold,10,
                return_all=True)

        return H,ransac_data['inliers']

3拼接图像：

估计出图像间的单应性矩阵(使用 RANSAC 算法)，现在我们需要将所有的图像扭曲到一个公共的图像平面上。通常，这里的公共平面为中心图像平面(否则，需要进行大量变形)。一种方法是创建一个很大的图像，比如图像中全部填充 0，使其和中心图像平行，然后将所有的图像扭曲到上面。由于我们所有的图像是由照相机水平旋转拍摄的，因此我们可以使用一个较简单的步骤:将中心图像左边或者右边的区域填充0，以便为扭曲的图像腾出空间。

def panorama(H,fromim,toim,padding=2400,delta=2400):
    """ 使用单应性矩阵 H(使用 RANSAC 健壮性估计得出)，协调两幅图像，创建水平全景图像。结果
        为一幅和 toim 具有相同高度的图像。padding 指定填充像素的数目，delta 指定额外的平移量
    """

    # 检查图像是灰度图像，还是彩色图像

    is_color = len(fromim.shape) == 3
    # 用于 geometric_transform() 的单应性变换

    def transf(p):
        p2 = dot(H,[p[0],p[1],1])
        return (p2[0]/p2[2],p2[1]/p2[2])

    if H[1,2]<0: # fromim 在右边
        print 'warp - right'
        # 变换 fromim
        if is_color:
            # 在目标图像的右边填充 0
            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:
            # 在目标图像的右边填充 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 = 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