计算机视觉-图像到图像的映射

席惜兮兮

已于 2022-05-25 22:03:40 修改

阅读量2.6k

点赞数

文章标签： python

于 2022-04-14 13:23:06 首次发布

本文链接：https://blog.csdn.net/Huanghuali_/article/details/124168297

版权

一、单应性变化

单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。在这里，平面是指图像或者三维中的平面表面。单应性变换具有很强的实用性，比如图像配准、图像纠正和纹理扭曲，以及创建全景图像。本质上，单应性变换 H，按照下面的方程映射二维中的点（齐次坐标意义下）：

或 $x{}'=Hx$

对于图像平面内的点（甚至是三维的），齐次坐标是个非常有用的表示方法。点的齐次坐标是依赖于其尺度定义的, x = [ x, y, w] = [αx, αy, αw] = [x/w,y/w,1] 都表示同一个二维点。因此，单应性矩阵H 也仅依赖尺度定义，所以，单应性矩阵具有 8 个独立的自由度。我们通常使用 w=1 来归一化点，这样，点具有唯一的图像坐标x和 y。这个额外的坐标使得可以简单地使用一个矩阵来表示变换。

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]))))

进行点和变换的处理时。我们会按照列优先的原则存储这些点。因此，n个二维点集将会存储为齐次坐标意义下的一个3×n数组。这种格式使得矩阵乘法和点的变换操作更加容易。对于其他的例子，比如对于聚类和分类的特征，我们将使用典型的行数组来存储数据。

在这些投影变换中，有一些特别重要的变换。比如，仿射变换：

或

保持了w=1，不具有投影变换所具有的强大变形能力，反射变换包括一个可逆矩阵A和一个平移向量 $t=\sqsubset t_x,t_y\sqsupset$ 。仿射变换可以用于很多应用，比如 图像扭曲。

相似变换：或

是一个包含尺度变化的二维刚体变换。上式中的向量s指定了变换的尺度，R是角度为θ的旋转矩阵， $t=\sqsubset t_x,t_y\sqsupset$ 在这里也是一个平移向量。如果s=1，那么该变换能够保持距离不变。此时，变换称为刚体变换。相似变换可以用于很多应用，比如图像配准。

下面让我们来一起探讨如何设计用于估计单应性矩阵的算法，然后看一下使用仿射变换进行图像扭曲，使用相似变换进行图像匹配，以及使用完全投影变换进行创建全景图像的一些例子。

1.直接线性变换算法

单应性矩阵可以有两幅图像（或者平面）中对应点对计算出来。前面已经提到过，一个完全射影变换具有8个自由度。根据对应点约束，每个对应点对可以写出两个方程，分别对应于x和y坐标。因此，计算单应性矩阵H需要4个对应点对。

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

或者Ah=0，其中A是一个具有对应点对二倍数量行数的矩阵。将这些对应点对方程的系数堆叠到一个矩阵红，我们可以使用SVD算法找到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]

上面的函数的第一步操作是检查点对两个数组中点的数目是否相同。如果不相同，函数将会抛出异常信息。这对于写出稳健的代码来说非常有用。

代码先对这些点进行归一化操作，使其均值为0，方差为1。因为算法的稳定性取决于坐标的表示情况和部分数值计算的问题，所以归一化操作非常重要。接下来我们使用对应点对来构造矩阵A。最小二乘解即为矩阵SVD分解后所得矩阵V的最后一行。该行经过变换后得到矩阵H。然后对这个矩阵进行处理和归一化，返回输出。

2.仿射变换

由于仿射变换具有6个自由度，因此我们需要三个对应点来估计矩阵H。通过将最后两个元素设置为0，即 $h_7=h_8=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著的Multiplr 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]

二、图像扭曲

对图像块应用仿射变换，我们将其称为图像扭曲（或者 仿射扭曲）。该操作不仅经常在计算机图形学中，而且经常出现在计算机视觉算法总。扭曲的操作可以使用SciPy工具包中的ndimage包来简单完成。命令：

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

使用上面所示的一个线性变换A和一个平移向量b来对图像块应用放射变换。选项参数 size 可以用来指定输出图像的大小。默认输出图像设置为和原始图像同样大小。我们可以运行下列代码来研究该函数是如何工作的

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

im = array(Image.open(r'guangqian.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]))

gray()
subplot(121)
imshow(im)
axis('off')
subplot(122)
imshow(im2)
axis('off')
show()

输出图像结果中丢失的像素用零补充

原始图像（左边）以及用 ndimage.affine_transform（）扭曲的图像

1.图像中的图像

仿射扭曲的一个简单例子是，将图像或者图像的一部分放置在另一幅图像中，是的他们能够和指定的区域或者标记物对齐。

将函数 image_in_image() 添加到 wary.py 文件中。该函数的输入参数为两幅图像和一个坐标。该坐标为将第一幅图像放置到第二幅图像中的角点坐标：

#图像中的图像
# -*- coding: utf-8 -*-
# import matplotlib.delaunay as md
from scipy.spatial import Delaunay
from scipy import ndimage
from pylab import *
from numpy import *
from PIL import Image
import homography

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]])
    # n是im2投影到im1的图像的长，m是宽
    # im1四个顶点的坐标

    # 得到单应性矩阵
    H = homography.Haffine_from_points(tp, fp)
    # 两个图形的变换关系（位移+缩放），只取H的前两行
    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

将扭曲的图像和第二幅图像融合，就创建 alpha 图像。该图像定义了每个像素从各个图像中获取的像素值成分多少。这里基于以下事实：扭曲的图像是在扭曲区域边界之外以 0 来填充的图像，来创建一个二值的 alpha 图像。严格意义上，需要在第一幅图象中的潜在 0 像素上加上一个小的数值，或者合理的处理这些 0 像素。这里我们使用图像坐标是齐次坐标意义下的。

使用仿射变换将一个表情包图像放置到另一幅图像中。

import warp
from numpy import *
from matplotlib.pyplot import *
from PIL import Image
import homography

# 仿射扭曲im1到im2的例子
im1 = array(Image.open(r'sese.jpg').convert('L'))
im2 = array(Image.open(r'guangqian.jpg').convert('L'))

gray()
subplot(131)
imshow(im1)
axis('equal')
axis('off')
subplot(132)
imshow(im2)
axis('equal')
axis('off')

# 选定4个点的y，x坐标
tp = array([[264,538,540,264],[40,36,605,605],[1,1,1,1]])

im3 = warp.image_in_image(im1, im2, tp)
subplot(133)
imshow(im3)
axis('equal')
axis('off')
show()

将两个图像分成两个三角形，然后分别进行扭曲图像操作

from matplotlib.pyplot import *
from numpy import array
from PIL import Image
from scipy import ndimage
import warp
import homography

im1 = array(Image.open(r'guangqian.jpg').convert('L'))
im2 = array(Image.open(r'shangda.jpg').convert('L'))

# 选定 im1 角上的一些点
m,n = im1.shape[:2]
fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])

# 选定一些目标点
tp = array([[50, 188, 175, 33], [60, 50, 375, 370], [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])
im3 = (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)*im3 + alpha*im1_t

figure()
gray()
imshow(im4)
axis('equal')
axis('off')
show()

2.分段仿射扭曲

三角形图像块的仿射扭曲可以完成角点的精确匹配。让我们看一下对应点对集合之间最常用的扭曲方式：分段仿射扭曲。给定任意图像的标记点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，我们可以将图像和另一幅图像的对应标记点扭曲对应。对于任何图形和图像处理库来说，这些都是最基本的操作。
为了三角化这些点，我们经常使用狄洛克三角剖分方法。

from numpy import *
from matplotlib.pyplot import *
from scipy.spatial import Delaunay

x, y = array(random.standard_normal((2, 100)))
tri = Delaunay(np.c_[x, y]).simplices

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')

figure()
plot(x, y, '*')
axis('off')
show()

随机二维点集的狄洛克三角部分示例

以下是用分段仿射图像扭曲的通用扭曲函数

# 打开图片，写入txt文件
from PIL import Image
from pylab import *
from warp import *

# 解决中文乱码
rcParams['font.sans-serif'] = 'SimHei'
rcParams['axes.unicode_minus'] = False

# 读取图片
im = array(Image.open('shangda.jpg'))

title('请按从左到右，从上到下的顺序点击30次：')
imshow(im)

points = ginput(30)

# 转整数
int_points = np.int_(points)

show()

# 写入txt文件
with open('points.txt', 'w') as f:
    for i in range(30):
        f.write(str(int_points[i][0]))
        f.write(" ")
        f.write(str(int_points[i][1]))
        f.write("\n")



 # 读取txt
 # -*- coding: utf-8 -*-
from pylab import *
from PIL import Image
import warp

# 打开图像，并将其扭曲
fromim = array(Image.open('shangda.jpg'))
x, y = meshgrid(range(5), range(6))
x = (fromim.shape[1]/4) * x.flatten()
y = (fromim.shape[0]/5) * y.flatten()

# 三角剖分
tri = warp.triangulate_points(x, y)

# 打开图像和目标点
im = array(Image.open('guangqian.jpg'))
tp = loadtxt('points.txt') # destination points

figure()
subplot(1, 4, 1)
axis('off')
imshow(im)

# 将点转换成齐次坐标
fp = array(vstack((y, x, ones((1, len(x))))), 'int')
tp = array(vstack((tp[:, 1], tp[:, 0], ones((1, len(tp))))), 'int')

# 扭曲三角形
im = warp.pw_affine(fromim, im, fp, tp, tri)

# 绘制图像
subplot(1, 4, 2)
axis('off')
imshow(fromim)
warp.plot_mesh(fp[1], fp[0], tri)
subplot(1, 4, 3)
axis('off')
imshow(im)
subplot(1, 4, 4)
axis('off')
imshow(im)
warp.plot_mesh(tp[1], tp[0], tri)
show()

右图（第一幅为带有标记物的目标图像，第二幅为带有三角剖分的图像，第三幅为扭曲后的图像，第四幅为带有三角剖分的扭曲图像）

三、创建全景图

此处记得把vl.dll ， vl.lib和sift.exe放到项目文件里面，记得改sift.py里面的cmmd路径

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

对于通用的geometric_transform()函数，我们需要指定能够描述像素到像素间映射的函数。在这个例子中，transf()函数就是该指定函数。干函数通过将像素和H相乘，然后对齐次坐标进行归一化来实现像素间的映射。通过查看H中的平移量，我们可以决定应该将图像填补到左边还是右边。当图像填补到左边时，由于目标图像中点的坐标也变化了，所以在“左边“情况中，需要在单应性矩阵中加入平移。简单起见，我们同样使用0像素的技巧来寻找alpha图。

1.拼接图像


# -*- coding: utf-8 -*-
from pylab import *
from numpy import *
from PIL import Image

# If you have PCV installed, these imports should work
import homography, warp
import sift

np.seterr(invalid='ignore')
"""
This is the panorama example from section 3.3.
"""

# set paths to data folder
featname = ['D:/PyCharmProjects/img/' + str(i + 1) + '.sift' for i in range(5)]  # 图片路径记得修改
imname = ['D:/PyCharmProjects/img/' + str(i + 1) + '.jpg' for i in range(5)]

# extract features and match
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])

# visualize the matches (Figure 3-11 in the book)
for i in range(4):
    im1 = array(Image.open(imname[i]))
    im2 = array(Image.open(imname[i + 1]))
    figure()
    sift.plot_matches(im2, im1, l[i + 1], l[i], matches[i], show_below=True)


# function to convert the matches to hom. points
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)

    # switch x and y - TODO this should move elsewhere
    fp = vstack([fp[1], fp[0], fp[2]])
    tp = vstack([tp[1], tp[0], tp[2]])
    return fp, tp


# estimate the homographies
model = homography.RansacModel()

fp, tp = convert_points(1)
H_12 = homography.H_from_ransac(fp, tp, model)[0]  # im 1 to 2

fp, tp = convert_points(0)
H_01 = homography.H_from_ransac(fp, tp, model)[0]  # im 0 to 1

tp, fp = convert_points(2)  # NB: reverse order
H_32 = homography.H_from_ransac(fp, tp, model)[0]  # im 3 to 2

tp, fp = convert_points(3)  # NB: reverse order
H_43 = homography.H_from_ransac(fp, tp, model)[0]  # im 4 to 3

# warp the images
delta = 1600  # for padding and translation

im1 = array(Image.open(imname[1]), "uint8")
im2 = array(Image.open(imname[2]), "uint8")
im_12 = warp.panorama(H_12, im1, im2, delta, delta)

im1 = array(Image.open(imname[0]), "f")
im_02 = warp.panorama(dot(H_12, H_01), im1, im_12, delta, delta)

im1 = array(Image.open(imname[3]), "f")
im_32 = warp.panorama(H_32, im1, im_02, delta, delta)

im1 = array(Image.open(imname[4]), "f")
im_42 = warp.panorama(dot(H_32, H_43), im1, im_32, delta, 2 * delta)

figure()
imshow(array(im_02, "uint8"))
axis('off')
savefig("HHL.png", dpi=300)
show()

为什么会出现黑色的呢，我用了几组图，都是一样的效果。要是有大佬知道，望赐教（抱拳）

（后续）啊我终于知道是什么问题了！！！！是因为排图的顺序。从右到左依次为1~5。

这是新的测试结果

成功了！！！（虽然还是有黑色的部分）

出现黑色可能是因为我把一张图片裁剪成五份进行拼接，或许做实验需要完全独立的图片

席惜兮兮

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
1
评论
计算机视觉-图像到图像的映射

一、单应性变化单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。在这里，平面是指图像或者三维中的平面表面。单应性变换具有很强的实用性，比如图像配准、图像纠正和纹理扭曲，以及创建全景图像。本质上，单应性变换 H，按照下面的方程映射二维中的点（齐次坐标意义下）：或对于图像平面内的点，齐次坐标是个非常有用的表示方法。点的齐次坐标是依赖于其尺度定义的, x = [ x, y, w] = [αx, αy, αw] = [x/w,y/w,1] 都表示同一个二维点。因此，单应性矩阵H ...
复制链接

扫一扫