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

最新推荐文章于 2024-08-31 21:39:38 发布

adchloe

最新推荐文章于 2024-08-31 21:39:38 发布

阅读量1.1k

点赞数 12

文章标签： python 计算机视觉算法

本文链接：https://blog.csdn.net/adchloe/article/details/141599539

版权

1 单应性变换

单应性变换是将一个平面内的点映射到另一个平面内对的二维投影变换。该变换可以用于图像匹配，图像纠正和纹理扭曲，以及创建全景图像。

单应性变换H，按照下面的方程映射二维中的点：
$\begin{bmatrix}x'\\y'\\w'\end{bmatrix}=\begin{bmatrix}h_1&h_2&h_3\\h_4&h_5&h_6\\h_7&h_8&h_9\end{bmatrix}\begin{bmatrix}x\\y\\w\end{bmatrix}\quad\text{或}\quad x'=Hx$

点的齐次坐标是依赖于其尺度定义的， $\mathbf{x}=[x,y,w]=[\alpha x,\alpha y,\alpha w]=[x/w,y/w,1]$ 都表示同一个二维点。单应性矩阵H仅依赖尺度定义，且具有8个独立的自由度。可以使用w=1来归一化点。

可以使用下面的函数来实现对点进行归一化和转换齐次坐标。

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.仿射变换
$\begin{bmatrix}x'\\y'\\1\end{bmatrix}=\begin{bmatrix}a_1&a_2&t_x\\a_3&a_4&t_y\\0&0&1\end{bmatrix}\begin{bmatrix}x\\y\\1\end{bmatrix}\quad\text{或}\quad x'=\begin{bmatrix}A&t\\0&1\end{bmatrix}x$
不具有投影变换所具有的强大变形能力，可用于图像扭曲
2.相似变换
$\begin{bmatrix}x'\\y'\\1\end{bmatrix}=\begin{bmatrix}s\cos(\theta)&-s\sin(\theta)&t_x\\s\sin(\theta)&s\cos(\theta)&t_y\\0&0&1\end{bmatrix}\begin{bmatrix}x\\y\\1\end{bmatrix}\quad\text{或}\quad x'=\begin{bmatrix}s\boldsymbol{R}&t\\\boldsymbol{0}&1\end{bmatrix}\boldsymbol{x}$
s=1时，变换能够保持距离不变，变换为刚体变换，可以用于图像配准。

1.1 直接线性变换算法

单应性矩阵可以由两幅图像中对应点对计算出来，由于完全射影变换具有8个自由度，计算单应性矩阵H徐娅4个对应点对。
DLT是给定4个或者更多对应点对矩阵，来计算单应性矩阵H的算法，可得到下面的方程：
$\begin{bmatrix}-x_1&-y_1&-1&0&0&0&x_1x'_1&y_1x'_1&x'_1\\0&0&0&-x_1&-y_1&-1&x_1y'_1&y_1y'_1&y'_1\\-x_2&-y_2&-1&0&0&0&x_2x'_2&y_2x'_2&x'_2\\0&0&0&-x_2&-y_2&-1&x_2y'_2&y_2y'_2&y'_2\\&\vdots&&\vdots&\vdots&&\vdots&\vdots\end{bmatrix}\begin{bmatrix}h_1\\h_2\\h_3\\h_4\\h_5\\h_6\\h_7\\h_8\\h_9\end{bmatrix}=\mathbf{0}$

或者Ah=0，可以使用SVD算法找到H的最小二乘解。代码函数如下：

def H_from_points(fp, tp):
    if fp.shape != tp.shape:
        raise RuntimeError("number of points don't match")
    m = np.mean(fp[:2],axis=1)
    maxstd = max(np.std(fp[:2],axis=1))+1e-9
    C1 = np.diag([1/maxstd,1/maxstd,1])
    C1[0][2] = -m[0]/maxstd
    C1[1][2] = -m[1]/maxstd
    fp = np.dot(C1,fp)
 
    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))
    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)
    H = V[8].reshape((3,3))
    H = np.dot(np.linalg.inv(C2),np.dot(H,C1))
    return H / H[2,2]

对这些点进行归一化操作，使其均值为0，方差为1。

1.2 仿射变换

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

下面函数使用对应点对来计算仿射变换矩阵：

def Haffine_from_points(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_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)

  A = concatenate((fp_cond[:2],tp_cond[:2]), axis=0)
  U,S,V = linalg.svd(A.T)
  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]

该算法类似于DLT算法，这些点需要经过预处理和去处理话操作。

2 图像扭曲

对图像块应用仿射变换，称为图像扭曲。可以使用Scipy工具包中的ndimage包来简单完成。
可以使用下述命令来对图像块应用仿射变换

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

代码如下：

from PIL import Image
from matplotlib import pyplot as plt
import numpy as np
from scipy import ndimage
 
img = np.array(Image.open(r'D:\test\pil.png').convert('L'))
H = np.array([[1.4,0.05,-100],[0.05,1.5,-100],[0,0,1]])
img2 = ndimage.affine_transform(img,H[:2,:2],(H[0,2],H[1,2]))
 
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.subplot(121), plt.imshow(img, plt.cm.gray), plt.title('原始图像')
plt.subplot(122), plt.imshow(img2, plt.cm.gray), plt.title('扭曲图像')
plt.show()

在这里插入图片描述
右图扭曲图像中丢失的像素使用0来填充。

2.1 图像中的图像

仿射扭曲的一个例子，将图像或者图像的一部分放置在另一幅图像中，使得它们能够和指定的区域或者标记物对齐。
代码如下：

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


im1 = array(Image.open('jimei.jpg').convert('L'))
im2 = array(Image.open('jimei2.jpg').convert('L'))

tp = array([[300,1250,1250,300],[0,0,280,280],[1,1,1,1]])
im3 = image_in_image(im1,im2,tp)

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.subplot(131), plt.imshow(im1, plt.cm.gray), plt.title('原始图像1'), plt.axis('off')
plt.subplot(132), plt.imshow(im2, plt.cm.gray), plt.title('原始图像2'), plt.axis('off')
plt.subplot(133), plt.imshow(im3, plt.cm.gray), plt.title('生成图像'), plt.axis('off')
plt.show()

将扭曲的图像和第二幅图像融合，创建了alpha图像，该图像定义了每个像素从各个图像中获取的像素值成分多少。

对于三个点，仿射变换可以将一幅图像进行扭曲，使这三对对应点对可以完美匹配。

2.2分段仿射扭曲

三角形图像块的仿射扭曲可以完成角点的精确匹配。给定任意图像的标记点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，可以将图像和另一幅图像的对应标记点扭曲对应。
使用狄洛克三角剖分方法，代码如下：

import numpy as np  
from scipy.spatial import Delaunay  
import matplotlib.pyplot as plt  
import random


def create_data():
    n = 50
    x = [random.randint(1, 50) for _ in range(n)]
    y = [random.randint(51, 100) for _ in range(n)]
    points = [(i, j) for i, j in zip(x, y)]
    points = np.array(points)
    return points


def create_delauney(points):
    # create a Delauney object using (x, y)
    tri = Delaunay(points)

    # paint a triangle
    plt.triplot(points[:, 0], points[:, 1], tri.simplices.copy(), c='black')
    plt.plot(points[:, 0], points[:, 1], 'o', c='green')
    plt.axis('equal')
    plt.show()


def main():
    # create random data
    point = create_data()
    # create your Delauney using random data
    create_delauney(point)


if __name__ == '__main__':
    main()

在这里插入图片描述
狄洛克剖分选择一些三角形，使三角剖分中所有三角形的最小角度最大。

2.3 图像配准

图像配准是对图像进行变换，使变换后的图像能够在常见的坐标系中对齐。配准分为严格配准和非严格配准。

3 创建全景图

在同一位置拍摄的两幅或多幅图像是单应性相关的。经常使用约束将很多图像缝补起来，拼称一个大的图像来创建全景图像。

3.1 RANSAC（随机一致性采样）

该方法用来找到正确模型来你和带有噪声数据的迭代方法。基本思想：数据中包含正确的点和噪声点，合理的模型能够在描述正确数据点的同时摒弃噪声点。

import numpy
import scipy  unavailable
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]  
        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:
   

    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)  
        return err_per_point


def test():
   
    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))  
    B_exact = scipy.dot(A_exact, perfect_fit)
    assert B_exact.shape == (n_samples, n_outputs)

   
    A_noisy = A_exact + numpy.random.normal(size=A_exact.shape)
    B_noisy = B_exact + numpy.random.normal(size=B_exact.shape)

    if 1:
        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))

   
    all_data = numpy.hstack((A_noisy, B_noisy))
    input_columns = range(n_inputs)  
    output_columns = [n_inputs + i for i in range(n_outputs)] 
    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])

  
    ransac_fit, ransac_data = ransac(all_data, model,
                                     5, 5000, 7e4, 50,  
                                     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()

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

SIFT具有很强的稳健性的描述子，能够比其他描述子产生更少的错误的匹配。

3.3 拼接图像

在估计出图像间的单应性矩阵后，需要将所有的图像扭曲到一个公共的图像平面上。
思想：创建一个很大的图像，使其和中心图像平行，然后将所有的图像扭曲到上面。

adchloe

关注

12
点赞
踩
17

收藏

觉得还不错? 一键收藏
0
评论
Python计算机视觉编程——第三章图像到图像的映射

给定任意图像的标记点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，可以将图像和另一幅图像的对应标记点扭曲对应。将扭曲的图像和第二幅图像融合，创建了alpha图像，该图像定义了每个像素从各个图像中获取的像素值成分多少。仿射扭曲的一个例子，将图像或者图像的一部分放置在另一幅图像中，使得它们能够和指定的区域或者标记物对齐。可以使用w=1来归一化点。在估计出图像间的单应性矩阵后，需要将所有的图像扭曲到一个公共的图像平面上。思想：创建一个很大的图像，使其和中心图像平行，然后将所有的图像扭曲到上面。
复制链接

扫一扫