Python将矩阵变换形状 python将矩阵转为图像

 

目录

引言

一、 单应性变换

1.1 直接线性变换算法

1.2 仿射变换 

二、 图像扭曲 

2.1 图像中的图像

2.2 分段仿射扭曲

2.2 图像配准

三、创建全景图

3.1 RANSAC（随机一致性采样）

3.2 拼接图像

引言

        本章讨论图像之间的变换，以及一些计算变换的方法，这些变换可以用于图像扭曲变形和图像配准。

一、 单应性变换

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

对于图像平面内的点，齐次坐标是一个非常有用的表达方式。点的齐次坐标是依赖于其尺度定义的，为此

，表示的是同一个二维点。单应性矩阵仅依赖尺度定义，具有8个独立的自由度。创建一个能够实现对点进行归一化和转换齐次坐标的功能：

import numpy as np

def normalize_homogeneous(points):
    """
    对点进行归一化和转换齐次坐标

    Args:
        points: ndarray, shape (N, D) or (N, D+1)
            待归一化的点坐标，每行表示一个点，每列表示维度

    Returns:
        points_norm: ndarray, shape (N, D+1)
            归一化后的点，转换为齐次坐标表示
    """
    # 将 points 转换为浮点型
    points = points.astype(np.float64)
    
    # 确定点的维度
    is_homogeneous = (points.shape[1] == 3)
    if is_homogeneous:
        points = points[:, :-1] / points[:, -1:]
    N, D = points.shape
    
    # 计算均值和方差
    mean = np.mean(points, axis=0)
    var = np.var(points, axis=0)
    std = np.sqrt(var)
    
    # 构造归一化矩阵
    scale = np.diag(1 / std)
    translate = -mean / std
    transform = np.eye(D)
    transform[:D-1, :D-1] = scale
    transform[:D-1, -1] = translate
    
    # 应用归一化矩阵并转换为齐次坐标表示
    points_norm = np.ones((N, D+1))
    points_norm[:, :-1] = points
    points_norm = points_norm @ transform.T
    
    if is_homogeneous:
        return points_norm
    else:
        return np.hstack((points_norm, np.ones((N, 1))))

进行点和变换的处理时，按照列优先的原则存储这些点，故n个三维点集将会依次存储为齐次坐标意义下的一个3*n数组。在这些投影变换中存在一些特别重要的变换如：

仿射变换：

w=1，不具有投影变换所具有的强大变形能力。仿射变换包含一个可逆矩阵A和一个平移向量

。

相似变换：

包含尺度变换的二维刚体变换，上式的s向量指定了变换的尺度，R是角度为的旋转矩阵，是一个平移向量，s=1时，该变换甭管保存距离不变，此时变换为刚体变换。

1.1 直接线性变换算法

        单应性矩阵可以由两幅图像中对于点计算出来，我们已经了解到一个完全射影变换具有8个自由度，而每个对应点可以写出两个方程，对应x和y坐标，所以计算单应性矩阵H要4个对应点对。计算H的方程如下：

python实现的算法如下：

from numpy import *
from numpy.linalg import svd, inv

def H_from_points(fp, tp):
    """
    使用DLT方法，计算单应性矩阵H，使fp映射到tp，点自动进行归一化
    :param fp:
    :param tp:
    :return:
    """
    if fp.shape != tp.shape:
        raise RuntimeError('number of points do not match')
    
    # 归一化点坐标
    def normalize_points(pts):
        mean = pts.mean(axis=1).reshape(-1, 1)
        std = pts.std(axis=1).mean()
        T = diag([1/std, 1/std, 1])
        T[:2, 2] = -mean[:2, 0] / std
        return T @ pts
    
    fp_norm = normalize_points(fp)
    tp_norm = normalize_points(tp)
    
    # 构造矩阵A
    nbr_correrpondences = fp.shape[1]
    A = zeros((2 * nbr_correrpondences, 9))
    for i in range(nbr_correrpondences):
        x, y = fp_norm[:, i]
        u, v = tp_norm[:, i]
        A[2 * i] = [-x, -y, -1, 0, 0, 0, u * x, u * y, u]
        A[2 * i + 1] = [0, 0, 0, -x, -y, -1, v * x, v * y, v]
    
    # 使用SVD求解最小二乘问题，得到单应性矩阵H
    _, _, V = svd(A)
    H = V[-1].reshape((3, 3))
    
    # 反归一化H并返回
    H = inv(normalize_points(tp)) @ H @ normalize_points(fp)
    return H / H[2, 2]

函数的第一步操作是检查点对两个数组中点的数目是否相同，如果不相同函数将会给出错误提示，对这些点进行归一化操作，使其均值为，方差为 ，因为算法的稳定性取决于坐标的表示情况和部分数值计算的问题，因此归一化非常重要，之后利用对应点对构造矩阵A，最小二乘解即为矩阵SVD分解后得到矩阵V的最后一行。该行经过变形后得到矩阵H，并对得到的矩阵进行处理和归一化，返回输出。

1.2 仿射变换 

        先前提到过仿射变换是一种重要的投影变换，其具有6个自由度，因此需要3个对应点估计矩阵H，将最后两个元素设置为0，就可以使用DLT算法估计得出，下面给出python的算法实现：
 

from numpy import *
from numpy.linalg import inv

def estimate_affine_matrix(src_points, dst_points):
    """
    使用DLT方法，计算仿射矩阵，使src_points映射到dst_points，需要至少3个对应点
    :param src_points:
    :param dst_points:
    :return:
    """
    if src_points.shape != dst_points.shape:
        raise RuntimeError('number of points do not match')
    if src_points.shape[1] < 3:
        raise RuntimeError('at least 3 points are needed')
    
    # 构造矩阵A
    nbr_correspondences = src_points.shape[1]
    A = zeros((2 * nbr_correspondences, 6))
    for i in range(nbr_correspondences):
        x, y = src_points[:, i]
        u, v = dst_points[:, i]
        A[2 * i] = [x, y, 1, 0, 0, 0]
        A[2 * i + 1] = [0, 0, 0, x, y, 1]

    # 使用SVD求解最小二乘问题，得到仿射矩阵
    _, _, V = svd(A)
    H = V[-1].reshape((2, 3))
    
    # 将最后两个元素设置为0，得到仿射矩阵
    H[0, 2] = 0
    H[1, 2] = 0
    
    return H

 在上述代码中，我们首先检查输入的点数是否足够，然后根据输入点的坐标构造矩阵 A。接下来，我们使用 SVD 求解最小二乘问题，得到一个包含 6 个元素的向量，然后将其重塑为一个 2x3 的仿射矩阵。最后，我们将矩阵的最后两个元素设置为零，以获得仿射变换的矩阵形式。

二、 图像扭曲 

        对图像应用仿射变换就将其称为图像扭曲，扭曲操作可使用Scipy工具包里的ndimage包实现，具体命令为：

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

     实验：

从实验结果中可以看出，经过ndimage.affine_transform()函数方法处理过的图像发生了偏折，且丢失的像素使用0进行填充。

2.1 图像中的图像

        将一幅图像或图像的一部分放置在另一幅图像中，使得它们能够和指定的标记物对齐，这就是仿射扭曲的简单例子。

        下面就使用仿射变换将一幅图像放置到另一幅图像做出实验：

 实验时无法读取图片，原先的代码：

img1 = cv2.imread(r"C:\Users\psen2\Desktop\ku\集大.jpg", 0)
img2 = cv2.imread(r"C:\Users\psen2\Desktop\ku\集大2.jpg", 0)

 尝试使用其他的图像读取函数来读取图片，例如Pillow库的Image.open()函数，例如：

img1 = np.array(Image.open(r"C:\Users\psen2\Desktop\ku\集大.jpg").convert('L'))
img2 = np.array(Image.open(r"C:\Users\psen2\Desktop\ku\集大2.jpg").convert('L'))

注意，在使用Pillow库的Image.open()函数时，需要将返回的PIL图像对象转为NumPy数组，同时还需要使用.convert('L')函数将图像转为灰度图像。

2.2 分段仿射扭曲

        三角形图像块的仿射扭曲可以实现角点的精确匹配，分段仿射扭曲是对应点对集合之间常见的扭曲方式。通过给定任意图像的标记点，并将这些点进行三角剖分，然后实验仿射扭曲来扭曲每个三角形，我们可以将图像和另一幅图像的对应标记点扭曲对应。

        三角化点的操作一般实验狄洛克三角剖分方法，其实验代码如下：

在Matplotlib库中就含有狄洛克三角剖分，需要我们注意的是不是在Pylab库中，所以在导入包时应该为：

from matplotlib.pyplot import *

输出结果为：

对于一幅彩色图像，当我们需要对其进行图像扭曲时就需要额外进行考虑，下面的函数提供了一个方法： 

def pw_affine(fromim,toim,fp,tp,tri):
    """
    从图像中扭曲矩形图像块
    :param fromim: 将要扭曲的图像
    :param toim: 目标图像
    :param fp: 齐次坐标下，扭曲前的点
    :param tp: 齐次坐标下，扭曲后的点
    :param tri: 三角部分
    :return: 
    """
                
    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

分析：这里首先检查图像是灰度图像还是彩色图像，当为彩色图像时，就要分别对每个颜色通道进行扭曲处理， 我们知道对于每个三角形而言，仿射变换是唯一确定的，因此使用了homography包中的Haffine_from_points()函数对其进行处理。

2.2 图像配准

             图像配准就是将不同时间、成像设备或不同条件下（比如气候、照度、摄像位置和角度等）获取的两幅或多幅图像进行匹配、叠加的过程，它已经被广泛地应用于遥感数据分析、计算机视觉和图像处理。

       由于配准图像的多样性和各种类型的退化，不能设计出适合所有配准任务的通用方法。每种方法不仅要考虑图像之间假定的几何变形类型，还要考虑辐射变形和噪声损坏，所需配准的准确率和应用数据特征。一般图像配准的流程为：

        1、对两幅图像进行特征提取得到特征点；

        2、通过进行相似性度量找到匹配的特征点对；

        3、通过匹配的特征点对得到图像空间坐标变换参数；

        4、最后由坐标变换参数进行图像配准。

在配准的过程中特征提取尤为关键，准确的特征提取为特征匹配的成功进行提供了保障。因此，寻求具有良好不变性和准确性的特征提取方法，对于匹配精度至关重要。先前我们学习到的SIRT特征提取就不失为一个很好的方法。

三、创建全景图

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

3.1 RANSAC（随机一致性采样）

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

RANSAC的基本假设是： 
（1）数据由“局内点”组成，例如：数据的分布可以用一些模型参数来解释； 
（2）“局外点”是不能适应该模型的数据； 
（3）除此之外的数据属于噪声。 
局外点产生的原因有：噪声的极值；错误的测量方法；对数据的错误假设。 
RANSAC也做了以下假设：给定一组（通常很小的）局内点，存在一个可以估计模型参数的过程；而该模型能够解释或者适用于局内点。在估计参数模型时，通常使用p表示一些迭代过程中从数据集内随机选取出来的点均为局内点的概率；此时结果模型就很可能会有作用，所以p也表征了算法所产生有用结果的概率。

假设每个点是真正内群的概率为 w：

w = 内群的数目/（内群数目+外群数目）

通常我们不知道 w 是多少,

是所选择的n个点都是内群的机率, 

是所选择的n个点至少有一个不是内群的机率,

是表示重复 k 次都没有全部的n个点都是内群的机率, 假设算法跑 k 次以后成功的机率是p，那么 

我们可以通过P反算得到抽取次数K：

所以如果希望成功机率高：

        1、当n不变时，k越大，则p越大；

        2、当w不变时，n越大，所需的k就越大。

        3、通常w未知，所以n选小一点比较好。 

下面给出实验示例：

import numpy as np
import scipy.linalg as sl
import scipy as sp


def ransac(symbol, model, n, k, t, d, debug=False, return_all=False):
    """
    :param symbol: 样本点
    :param model: 假设的模型
    :param n: 生成模型所需的最小样本点
    :param k: 最大迭代次数
    :param t: 阈值
    :param d: 拟合较好时，需要的样本点最少的个数，当做阈值看待
    :param debug:
    :param return_all:
    :return:
    """
    iterations = 0
    bestfit = None
    besterr = np.inf  # 设置默认值
    best_inlier_idxs = None
    while iterations < k:
        maybe_idxs, test_idxs = random_partition(n, symbol.shape[0])
        # 获取maybe{_idxs行数据（Xi，Yi）
        maybe_inliers = symbol[maybe_idxs, :]
        # 若干行（xi,yi）数据点
        test_points = symbol[test_idxs]
        # 拟合模型
        maybemodel = model.fit(maybe_inliers)
        # 计算误差
        test_err = model.get_error(test_points, maybemodel)
        also_idxs = test_idxs[test_err < t]
        also_inliers = symbol[also_idxs, :]
        if debug:
            print('test_err.min()', test_err.min())
            print('test_err.max*(', test_err.max())
            print('np.mean(test_err)', np.mean(test_err))
            print('iteration %d:len(also_inliers) = %d' % (iterations, len(also_inliers)))
        if len(also_inliers > d):
            bettersymbol = np.concatenate((maybe_inliers, also_inliers))
            bettermodel = model.fit(bettersymbol)
            better_err = model.get_error(bettersymbol, bettermodel)
            # 将平均误差作为新误差
            this_err = np.mean(better_err)
            if this_err < besterr:
                bestfit = bettermodel
                besterr = this_err
                best_inlier_idxs = np.concatenate((maybe_idxs, also_idxs))
            iterations += 1
        if bestfit is None:
            raise ValueError("didn't met fit criteria")
        if return_all:
            return bestfit, {'inliers': best_inlier_idxs}
        else:
            return bestfit


def random_partition(n, n_symbol):
    # 获取n_symbol下的索引
    all_idxs = np.arange(n_symbol)
    # 打乱下标索引
    np.random.shuffle(all_idxs)
    idxs_1 = all_idxs[:n]
    idxs_2 = all_idxs[n:]
    return idxs_1, idxs_2


class LinearLeastSquaresModel:
    """
    最小二乘求线性解用于ransac的输入模型
    """

    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, symbol):
        # 第一列xi-->行xi
        A = np.vstack([symbol[:, i] for i in self.input_columns]).T
        # 第第二列yi-->行yi
        B = np.vstack([symbol[:, i] for i in self.output_columns]).T
        # residues：残差和
        x, resids, rank, s = sl.lstsq(A, B)
        # 返回最小平方和向量
        return x

    def get_error(self, symbol, model):
        A = np.vstack([symbol[:, i] for i in self.input_columns]).T
        B = np.vstack([symbol[:, i] for i in self.output_columns]).T
        # 计算y值B_fit = model.k+A+model.b
        B_fit = sp.dot(A, model)
        err_per_point = np.sum((B - B_fit) ** 2, axis=1)
        return err_per_point


def test():
    # 生成理想数据
    n_samples = 500
    n_inputs = 1
    n_output = 1
    A_exact = 20 * np.random.random((n_samples, n_inputs))
    perfect_fit = 60 * np.random.normal(size=(n_inputs, n_output))
    B_exact = sp.dot(A_exact, perfect_fit)
    # 加入高斯噪声
    A_noise = A_exact + np.random.normal(size=A_exact.shape)
    B_noise = B_exact + np.random.normal(size=B_exact.shape)

    if 1:
        # 添加局外点
        n_outliers = 100
        all_idxs = np.arange(A_noise.shape[0])
        np.random.shuffle(all_idxs)
        outlier_idxs = all_idxs[:n_outliers]
        non_outlier_idxs = all_idxs[n_outliers:]
        A_noise[outlier_idxs] = 20 * np.random.random((n_outliers, n_inputs))
        B_noise[outlier_idxs] = 50 * np.random.normal(size=(n_outliers, n_output))

    all_data = np.hstack((A_noise, B_noise))
    input_columns = range(n_inputs)  # the first columns of the array
    output_columns = [n_inputs + i for i in range(n_output)]  # the last columns of the array
    debug = False
    model = LinearLeastSquaresModel(input_columns, output_columns, debug=debug)

    linear_fit, resids, rank, s = np.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,
                                     debug=debug, return_all=True)
    if 1:
        import pylab

        sort_idxs = np.argsort(A_exact[:, 0])
        A_col0_sorted = A_exact[sort_idxs]

        if 1:
            pylab.plot(A_noise[:, 0], B_noise[:, 0], 'k.', label='data')
            pylab.plot(A_noise[ransac_data['inliers'], 0], B_noise[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],
                   np.dot(A_col0_sorted, ransac_fit)[:, 0],
                   label='RANSAC fit')
        pylab.plot(A_col0_sorted[:, 0],
                   np.dot(A_col0_sorted, perfect_fit)[:, 0],
                   label='exact system')
        pylab.plot(A_col0_sorted[:, 0],
                   np.dot(A_col0_sorted, linear_fit)[:, 0],
                   label='linear fit')
        pylab.legend()
        pylab.show()


if __name__ == '__main__':
    test()

 输出结果为：

分析：ransac就是用来解决样本中的噪声问题的，且最多可处理50%的噪声情况。其优点在于能够很好的估计模型参数，在运行过程中，如果出现了一种足够好的模型，就会跳出主循环，当直接从maybe_model计算this_err时虽然会节约比较两种模型的错误时间，但对于噪声会更加的敏感。不过由于其计算参数的迭代次数没有上限，所以其运行的时间没有保证，倘若强行设置迭代次数可能会导致错误的结果。

3.2 拼接图像

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

代码：

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取两张测试图像
img1 = cv2.imread('2.jpg')
img2 = cv2.imread('4.jpg')

# 显示输入的两张图像
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 6))
ax1.imshow(cv2.cvtColor(img1, cv2.COLOR_BGR2RGB))
ax1.set_axis_off()
ax1.set_title('Image 1')
ax2.imshow(cv2.cvtColor(img2, cv2.COLOR_BGR2RGB))
ax2.set_axis_off()
ax2.set_title('Image 2')
plt.tight_layout()
plt.show()

# 使用 SIFT 特征检测器和描述子提取器提取特征点和描述子
sift = cv2.SIFT_create()
kp1, des1 = sift.detectAndCompute(img1, None)
kp2, des2 = sift.detectAndCompute(img2, None)

# 使用暴力匹配器计算匹配关系
bf = cv2.BFMatcher()
matches = bf.match(des1, des2)

# 选取前100个匹配关系
matches = sorted(matches, key=lambda x: x.distance)[:100]

# 从匹配关系中提取对应的点坐标
x1 = np.array([kp1[m.queryIdx].pt for m in matches])
x2 = np.array([kp2[m.trainIdx].pt for m in matches])

# 使用 RANSAC 算法估计稳健的单应性矩阵
H, _ = cv2.findHomography(x2, x1, cv2.RANSAC, 5.0)

# 将图像2进行单应性变换，拼接到图像1上
h, w, _ = img1.shape
dst = cv2.warpPerspective(img2, H, (w, h))
print(f"dst shape: {dst.shape}")
print(f"img1 shape: {img1.shape}")
merged = cv2.addWeighted(img1, 0.5, dst, 0.5, 0)

# 显示拼接后的图像
fig, ax = plt.subplots(figsize=(12, 6))
ax.imshow(cv2.cvtColor(merged, cv2.COLOR_BGR2RGB))
ax.set_axis_off()
ax.set_title('Merged Image')
plt.tight_layout()
plt.show()

# 将输入的两张图像和拼接后的图像一起输出
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2, figsize=(12, 8))
ax1.imshow(cv2.cvtColor(img1, cv2.COLOR_BGR2RGB))
ax1.set_axis_off()
ax1.set_title('Image 1')
ax2.imshow(cv2.cvtColor(img2, cv2.COLOR_BGR2RGB))
ax2.set_axis_off()
ax2.set_title('Image 2')
ax3.imshow(cv2.cvtColor(merged, cv2.COLOR_BGR2RGB))
ax3.set_axis_off()
ax3.set_title('Merged Image')
ax4.set_axis_off()
plt.tight_layout()
plt.show()

输出结果如下：