[图像处理]图片的伸缩变换

要求

 编写程序，从标准输入读取图片的四个坐标和 bmp 格式图片，将这四个坐标框出的四边形拉伸成 224*224 的正方形，显示，并以 bmp 格式输出到标准输出。

思路

 给定图像的四个坐标 $\{x_i,y_i\} ,(i=1,2,3,4)$ 并且已知每个源坐标经变换后的对应坐标 $\{ x'_i,y'_i \} ,(i=1,2,3,4)$

普通线性变换

想到的第一个方法，对源坐标进行线性变换

$\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} x'\\ y' \end{bmatrix}$

但是很快会发现基本的线性变换并不能满足要求，因为普通的线性变换只能做到平行四边形到平行四边形的变换，并且不能应对图像的平移变换。

仿射变换

接着来看看仿射变换，仿射变换能够应对一般的图像旋转(线性变换)，平移(向量加)，缩放(线性变换)，错切，反转。
在二维下的仿射变换形式：

$\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} + \begin{bmatrix} b_{1} \\ b_{2} \end{bmatrix} = \begin{bmatrix} x'\\ y' \end{bmatrix}$

但是写成这样的形式并不利于求解。我们把上式展开写：

$\left\{\begin{matrix} a_{11}x + a_{12}y + b_1 = x'\\ a_{21}x + a_{22}y + b_2 = y' \end{matrix}\right.$

所以就可以先把二维的坐标换为三维齐次形式，并把仿射变换写成如下形式：

$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} = \begin{bmatrix} x' \\ y' \\ 0 \end{bmatrix}$

 具体如何求解还是回到展开式：$$\left{\begin{matrix} a_{11}x + a_{12}y + a_{13} = x’\ a_{21}x + a_{22}y + a_{23} = y’ \end{matrix}\right.$$ 发现这是一个六元方程组于是给定图像中三个不同源点坐标的对应变换后的坐标就可以求解仿射变换矩阵，我们把方程组写成如下形式：

$\left\{\begin{matrix} a_{11}x_1 + a_{12}y_1 + a_{13} + 0 + 0 + 0 = x'_1 \\ 0 + 0 + 0 + a_{21}x_1 + a_{22}y_1 + a_{23} = y'_1 \\ a_{11}x_2 + a_{12}y_2 + a_{13} + 0 + 0 + 0 = x'_2 \\ 0 + 0 + 0 + a_{21}x_2 + a_{22}y_2 + a_{23} = y'_2 \\ a_{11}x_3 + a_{12}y_3 + a_{13} + 0 + 0 + 0 = x'_3 \\ 0 + 0 + 0 + a_{21}x_3 + a_{22}y_3 + a_{23} = y'_3 \end{matrix}\right.$

用矩阵来表达即为：

$\begin{bmatrix} x_1 & y_1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & x_1 & y_1 & 1 \\ x_2 & y_2 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & x_2 & y_2 & 1 \\ x_2 & y_2 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & x_2 & y_2 & 1 \end{bmatrix} \begin{bmatrix} a_{11} \\ a_{12} \\ a_{13} \\ a_{21} \\ a_{22} \\ a_{23} \end{bmatrix} = \begin{bmatrix} x'_1 \\ y'_1 \\ x'_2 \\ y'_2 \\ x'_3 \\ y'_3 \end{bmatrix}$

记为：$XA=X'$ 。则,可以解出 $A = X^{-1}X'$

 但是此变换仍有局限性，因为本质上仿射变换只是在普通线性变换的基础上添加了平移量，所以原图像中平行四边形在变换后图像中仍是平行四边形。由于是平行四边形，只需要指定三个顶点即可确定。这显然还是不满足要求，寻找新的方法。

透视变换(投影变换)

 直观的理解透视变换，将二维平面中 $(X,Y)$ 转换为三维中 $(X,Y,Z)$，再到另外一个二维 $(X',Y')$ 的映射。相较于仿射变换，透视变换不仅能够做到平移，线性变换，还能够提供更多的自由度实现更为复杂的变换。

齐次化

 在进行透视变换之前，因为是将二维平面映射到三维中，仍需要像之前处理仿射变换一样将坐标齐次化。

$\begin{bmatrix} x \\ y \end{bmatrix} \rightarrow \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}$

变换方程

 透视变换的方程也就是将仿射变换进行扩展：

$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & 1 \end{bmatrix} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} = \begin{bmatrix} X \\ Y \\ Z \end{bmatrix}$

$\begin{bmatrix} x' \\ y' \\ 1 \end{bmatrix} = \begin{bmatrix} \frac{X}{Z} \\ \frac{Y}{Z} \\ 1 \end{bmatrix}$

 则可以将变换方程写为：

$\left\{ \begin{matrix} X = a_{11}x + a_{12}y + a_{13} \\ Y = a_{21}x + a_{22}y + a_{23} \\ Z = a_{31}x + a_{32}y + 1 \end{matrix} \right.$

$\left\{ \begin{matrix} x' = \frac{a_{11}x + a_{12}y + a_{13}}{a_{31}x + a_{32}y + 1} \\ y' = \frac{a_{21}x + a_{22}y + a_{23}}{a_{31}x + a_{32}y + 1} \\ 1 = \frac{Z}{Z} \end{matrix} \right.$

 我们将第二个方程组的两边同时乘以分母再移相，可以将方程组化成以下形式：

$\left\{ \begin{matrix} a_{11}x + a_{12}y + a_{13} - a_{31}xx' - a_{32}yx' = x' \\ a_{21}x + a_{22}y + a_{23} - a_{31}xy' - a_{32}yy' = y' \end{matrix} \right.$

因为$x,y,x',y'$都是已知量所以这是一个八元方程组，因此如果有4对坐标点我们就可以通过方程组解出投影变换矩阵。

$\left\{ \begin{matrix} a_{11}x_1 + a_{12}y_1 + a_{13} - a_{31}x_1x'_1 - a_{32}y_1x'_1 = x'_1 \\ a_{21}x_1 + a_{22}y_1 + a_{23} - a_{31}x_1y'_1 - a_{32}y_1y'_1 = y'_1 \\ a_{11}x_2 + a_{12}y_2 + a_{13} - a_{31}x_2x'_2 - a_{32}y_1x'_2 = x'_2 \\ a_{21}x_2 + a_{22}y_2 + a_{23} - a_{31}x_2y'_2 - a_{32}y_1y'_2 = y'_2 \\ a_{11}x_3 + a_{12}y_3 + a_{13} - a_{31}x_3x'_3 - a_{32}y_1x'_3 = x'_3 \\ a_{21}x_3 + a_{22}y_3 + a_{23} - a_{31}x_3y'_3 - a_{32}y_1y'_3 = y'_3 \\ a_{11}x_4 + a_{12}y_4 + a_{13} - a_{31}x_4x'_4 - a_{32}y_1x'_4 = x'_4 \\ a_{21}x_4 + a_{22}y_4 + a_{23} - a_{31}x_4y'_4 - a_{32}y_1y'_4 = y'_4 \end{matrix} \right.$

 这样就可以把方程组写成矩阵形式：

$\begin{bmatrix} x_1 & y_1 & 1 & 0 & 0 & 0 & -x_1x'_1 & -y_1x'_1 \\ 0 & 0 & 0 & x_1 & y_1 & 1 & -x_1y'_1 & -y_1y'_1 \\ x_2 & y_2 & 1 & 0 & 0 & 0 & -x_2x'_2 & -y_2x'_2 \\ 0 & 0 & 0 & x_2 & y_2 & 1 & -x_2y'_2 & -y_2y'_2 \\ x_3 & y_3 & 1 & 0 & 0 & 0 & -x_3x'_3 & -y_3x'_3 \\ 0 & 0 & 0 & x_3 & y_3 & 1 & -x_3y'_3 & -y_3y'_3 \\ x_4 & y_4 & 1 & 0 & 0 & 0 & -x_4x'_4 & -y_4x'_4 \\ 0 & 0 & 0 & x_4 & y_4 & 1 & -x_4y'_4 & -y_4y'_4 \end{bmatrix} \begin{bmatrix} a_{11} \\ a_{12} \\ a_{13} \\ a_{21} \\ a_{22} \\ a_{23} \\ a_{31} \\ a_{32} \end{bmatrix} = \begin{bmatrix} x'_1 \\ y'_1 \\ x'_2 \\ y'_2 \\ x'_3 \\ y'_4 \\ x'_4 \\ y'_4 \end{bmatrix}$

同样将上式记作$AB=X$ ，则可以解得 $B=A^{-1}X$

具体实现代码

编译依赖库

Eigen3(开源C++线性代数库),OpenCV3.8.4(其余版本未测试)。请在项目CMakeLists.txt文件中加入Eigen3和OpenCV依赖

代码

#include <iostream>
#include <string>
#include <vector>
#include <Eigen/Dense>

#include <opencv2/core/core.hpp>
#include <opencv2/features2d/features2d.hpp>
#include <opencv2/highgui/highgui.hpp>

using namespace std;
using namespace Eigen;
using namespace cv;

const int row_count = 224;
const int col_count = 224;

// 透射矩阵解算
Matrix3d getPerspectiveTransform(vector<KeyPoint> Src, vector<KeyPoint> Dst)
{
    if(Src.size() != Dst.size() && Src.size() != 4)
    {
        cerr << "input error" << endl;
        exit(1);
    }

    MatrixXd A(8, 8), B(8, 1);
    A.fill(0), B.fill(0); // 初始化矩阵

    // 初始化矩阵系数
    for(int i = 0; i < 4; ++i)
    {
        A(2*i, 0) = Src[i].pt.x, A(2*i, 1) = Src[i].pt.y, A(2*i, 2) = 1;
        A(2*i, 6) = -Src[i].pt.x*Dst[i].pt.x, A(2*i, 7) = -Src[i].pt.y*Dst[i].pt.x;

        A(2*i + 1, 3) = Src[i].pt.x, A(2*i + 1, 4) = Src[i].pt.y, A(2*i + 1, 5) = 1;
        A(2*i + 1, 6) = -Src[i].pt.x*Dst[i].pt.y, A(2*i + 1, 7) = -Src[i].pt.y*Dst[i].pt.y;

        B(2*i, 0) = Dst[i].pt.x, B(2*i + 1, 0) = Dst[i].pt.y;
    }

    MatrixXd result = A.inverse() * B;
    Matrix3d WarpMatrix;

    // 结果重新整形
    for(int i = 0; i < 3; i++)
        for(int j = 0; j < 3; j++)
        {
            if(i == 2 && j == 2) WarpMatrix(i, j) = 1.0;
            else WarpMatrix(i, j) = result(i * 3 + j, 0);
        }
    return WarpMatrix;
}

int main(int argc, char** argv) {

    string SrcFile;
    cout << "Please input the source image." << endl;
    cin >> SrcFile;

    Mat _img = imread(SrcFile, CV_LOAD_IMAGE_COLOR); // 读取图像
    Mat draw_img; // 目标点图像
    Mat out_img(row_count, col_count, CV_8UC3); // 输出图像

    vector<KeyPoint> Src, Dst; // 目标点

    // 点输入
    cout << "Please input points." << endl;
    for(int i = 0; i < 4; ++i)
    {
        int x, y;
        cin >> x >> y;
        Dst.emplace_back((double)x, (double)y, 10);
    }

    // 初始化目标坐标
    Src.emplace_back(0, 0, 10);Src.emplace_back(0, col_count - 1, 10);
    Src.emplace_back(row_count - 1, col_count - 1, 10);Src.emplace_back(row_count - 1, 0, 10);

    MatrixXd WarpMatrix = getPerspectiveTransform(Src, Dst); // 解算投影矩阵

    for(int i = 0; i < Dst.size(); ++i) // KeyPoint 的x y 坐标与实际坐标反转
        swap(Dst[i].pt.x, Dst[i].pt.y);

    drawKeypoints(_img, Dst, draw_img);

    // 计算投影坐标
    Vector3d s, t;
    s(2) = 1;
    for(int i = 0; i < row_count; ++i)
        for(int j = 0; j < col_count; ++j)
        {
            s(0) = i;
            s(1) = j;
            t = WarpMatrix * s;

            // 归一化
            t(0) /= t(2);
            t(1) /= t(2);

            out_img.at<Vec3b>(i, j) = _img.at<Vec3b>((uint)t(0), (uint)t(1));
        }

    imshow("Source Image", draw_img);
    imshow("Out Image", out_img);
    imwrite("OutImage.bmp", out_img);

    waitKey(0);
    return 0;
}