上篇 OpenCV 之圖象幾何變換介紹了等距、相似和仿射變換，本篇側重投影變換的平面單應性、OpenCV相關函式、應用實體等，

1 投影變換

1.1 平面單應性

投影變換 (Projective Transformation)，是仿射變換的泛化 (或普遍化)，二者區別如下：

假定平面 $P^{2}$ 與 $Q^{2}$ 之間，存在映射 $H_{3 \times 3}$，使得 $P^{2}$ 內任意點 $(x_p, y_q, 1)$，滿足下式：

$\quad \begin{bmatrix} x_q \\ y_q \\ 1 \end{bmatrix} = \begin{bmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{bmatrix} \begin{bmatrix} x_p \\ y_p \\ 1 \end{bmatrix} = H_{3 \times 3} \begin{bmatrix} x_p\\ y_p \\ 1 \end{bmatrix}$

當 $H$ 非奇異時，$P^{2}$ 與 $Q^{2}$ 之間的映射即為 2D 投影變換，也稱平面單應性， $H_{3 \times 3}$ 則為單應性矩陣

例如：在相機標定中，如果選取了二維平面標定板，則物平面和像平面之間的映射，就是一種典型的平面單應性

1.2 單應性矩陣

$H$ 有 9 個未知數，但實際只有 8 個自由度 (DoF)，其歸一化通常有兩種方法：

第一種方法，令 $h_{33}=1$；第二種方法，加單位向量限制 $h_{11}^2+h_{12}^2+h_{13}^2+h_{21}^2+h_{22}^2+h_{23}^2+h_{31}^2+h_{32}^2+h_{33}^2=1$

下面接著第一種方法，繼續推導公式如下：

$\quad x' = \dfrac{h_{11}x + h_{12}y + h_{13}}{h_{31}x + h_{32}y + 1} $

$\quad y' = \dfrac{h_{21}x + h_{22}y + h_{23}}{h_{31}x + h_{32}y + 1} $

整理得：

$\quad x \cdot h_{11} + y \cdot h_{12} + h_{13} - xx' \cdot h_{31} - yx' \cdot h_{32} = x' $

$\quad x \cdot h_{21} + y \cdot h_{22} + h_{23} - xy' \cdot h_{31} - yy' \cdot h_{32} = y' $

一組對應特征點 $(x, y) $ -> $ (x', y')$ 可構造 2 個方程，要求解 8 個未知數 (歸一化后的)，則需要 8 個方程，4 組對應特征點

$\quad \begin{bmatrix} x_{1} & y_{1} & 1 &0 &0&0 & -x_{1}x'_{1} & -y_{1}x'_{1} \\ 0&0&0& x_{1} & y_{1} &1& -x_{1}y'_{1} & -y_{1}y'_{1} \\ x_{2} & y_{2} & 1 &0 &0&0 & -x_{2}x'_{2} & -y_{2}x'_{2} \\ 0&0&0& x_{2} & y_{2} &1& -x_{2}y'_{2} & -y_{2}y'_{2} \\ x_{3} & y_{3} & 1 &0 &0&0 & -x_{3}x'_{3} & -y_{3}x'_{3} \\ 0&0&0& x_{3} & y_{3} &1& -x_{3}y'_{3} & -y_{3}y'_{3} \\ x_{4} & y_{4} & 1 &0 &0&0 & -x_{4}x'_{4} & -y_{4}x'_{4} \\ 0&0&0& x_{4} & y_{4} &1& -x_{4}y'_{4} & -y_{4}y'_{4} \end{bmatrix} \begin{bmatrix} h_{11} \\ h_{12}\\h_{13}\\h_{21}\\h_{22}\\h_{23}\\h_{31}\\h_{32} \end{bmatrix} = \begin{bmatrix} x'_{1}\\y'_{1}\\x'_{2}\\y'_{2}\\x'_{3}\\y'_{3}\\x'_{4}\\y'_{4} \end{bmatrix} $

因此，求 $H$ 可轉化為 $Ax=b$ 的通用解，參見 OpenCV 中 getPerspectiveTransform() 函式的原始碼實作

2 OpenCV 函式

2.1 投影變換矩陣

a) 四組對應特征點：已知四組對應特征點坐標，帶入 getPerspectiveTransform() 函式中，可求解 src 投影到 dst 的單應性矩陣 $H_{3 \times 3}$

Mat getPerspectiveTransform (
    const Point2f     src[],    // 原影像的四角頂點坐標
    const Point2f     dst[],    // 目標影像的四角頂點坐標
    int solveMethod = DECOMP_LU // solve() 的解法
)

該函式的源代碼實作如下：構造8組方程，轉化為 $Ax=b$ 的問題，呼叫 solve() 函式來求解

Mat getPerspectiveTransform(const Point2f src[], const Point2f dst[], int solveMethod)
{
    Mat M(3, 3, CV_64F), X(8, 1, CV_64F, M.ptr());
    double a[8][8], b[8];
    Mat A(8, 8, CV_64F, a), B(8, 1, CV_64F, b);

    for( int i = 0; i < 4; ++i )
    {
        a[i][0] = a[i+4][3] = src[i].x;
        a[i][1] = a[i+4][4] = src[i].y;
        a[i][2] = a[i+4][5] = 1;
        a[i][3] = a[i][4] = a[i][5] =
        a[i+4][0] = a[i+4][1] = a[i+4][2] = 0;
        a[i][6] = -src[i].x*dst[i].x;
        a[i][7] = -src[i].y*dst[i].x;
        a[i+4][6] = -src[i].x*dst[i].y;
        a[i+4][7] = -src[i].y*dst[i].y;
        b[i] = dst[i].x;
        b[i+4] = dst[i].y;
    }

    solve(A, B, X, solveMethod);
    M.ptr<double>()[8] = 1.;

    return M;
}

b) 多組對應特征點：對于兩個平面之間的投影變換，只要求得對應的兩組特征點，帶入 findHomography() 函式，便可得到 srcPoints 投影到 dstPoints 的 $H_{3 \times 3}$

Mat findHomography (
    InputArray      srcPoints,   // 原始平面特征點坐標，型別是 CV_32FC2 或 vector<Point2f> 
    InputArray      dstPoints,   // 目標平面特征點坐標，型別是 CV_32FC2 或 vector<Point2f> 
    int method = 0,              // 0--最小二乘法; RANSAC--基于ransac的方法
    double ransacReprojThreshold = 3, // 最大允許反投影誤差
    OutputArray mask = noArray(),     // 
    const int maxIters = 2000,        // 最多迭代次數
    const double confidence = 0.995   // 置信水平
)

2.2 投影變換影像

已知單應性矩陣$H_{3 \times 3}$，將任意影像代入 warpPerspective() 中，便可得到經過 2D投影變換的目標影像

void warpPerspective (
    InputArray      src,  // 輸入影像
    OutputArray     dst,  // 輸出圖象(大小 dsize，型別同 src)
    InputArray        M,  // 3x3 單應性矩陣
    Size          dsize,  // 輸出影像的大小
    int     flags = INTER_LINEAR,         // 插值方法 
    int     borderMode = BORDER_CONSTANT, //
    const Scalar& borderValue = https://www.cnblogs.com/xinxue/p/Scalar()  // 
)

3 單應性的應用

3.1 影像校正

像平面 image1 和 image2 是相機在不同位置對同一物平面所成的像，分別對應右下圖的 PP1、PP2 和 PP3，這三個平面任選兩個都互有平面單應性

以相機標定程序為例，當相機從不同角度對標定板進行拍攝時，利用任意兩個像平面之間的單應性，可將各角度拍攝的影像，轉換為某一特定視角的影像

3.2 影像拼接

當相機圍繞其投影軸，只做旋轉運動時，所有的像素點可等效視為在一個無窮遠的平面上，則單應性矩陣可由旋轉變換 $R$ 和相機標定矩陣 $K$ 來表示

$\quad \large s \begin{bmatrix} x' \\ y' \\1 \end{bmatrix} = \large K \cdot \large R \cdot \large K^{-1} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}$

因此，如果已知相機的標定矩陣，以及旋轉變換前后的位姿，可利用平面單應性，將旋轉變換前后的兩幅影像拼接起來

3.3 位姿估計

如下所示，當相機對著一個帶多個特征點的平面拍攝時，物、像平面之間便有了單應性映射 $H_1$：定義 $\widehat{n}$ 為平面的法向量，$d$ 為物平面到相機的距離 (沿著$\widehat{n}$的方向)

如果相機變換位姿，從不同角度對該平面進行成像，則可建立起相機所有位姿和該平面的單應性映射 $H_2，H_3，H_4 ... $ ，從而計算得出相機的位姿 (也稱 PnP 問題)

4 代碼示例

4.1 影像校正

#include <iostream>

#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/highgui.hpp"
#include "opencv2/calib3d.hpp"

using namespace std;
using namespace cv;

Size kPatternSize = Size(9, 6);

int main()
{
    // 1) read image
    Mat src = https://www.cnblogs.com/xinxue/p/imread("chessboard1.jpg");
    Mat dst = imread("chessboard2.jpg");
    if (src.empty() || dst.empty())
        return -1;

    // 2) find chessboard corners
    vector<Point2f> corners1, corners2;
    bool found1 = findChessboardCorners(src, kPatternSize, corners1);
    bool found2 = findChessboardCorners(dst, kPatternSize, corners2);
    if (!found1 || !found2)
        return -1;

    // 3) get subpixel accuracy of corners
    Mat src_gray, dst_gray;
    cvtColor(src, src_gray, COLOR_BGR2GRAY);
    cvtColor(dst, dst_gray, COLOR_BGR2GRAY);
    cornerSubPix(src_gray, corners1, Size(11, 11), Size(-1, -1), TermCriteria(TermCriteria::EPS | cv::TermCriteria::COUNT, 30, 0.1));
    cornerSubPix(dst_gray, corners2, Size(11, 11), Size(-1, -1), TermCriteria(TermCriteria::EPS | cv::TermCriteria::COUNT, 30, 0.1));

    // 4)
    Mat H = findHomography(corners1, corners2, RANSAC);
    // 5) 
    Mat src_warp;
    warpPerspective(src, src_warp, H, src.size());

    // 6) 
    imshow("src", src);
    imshow("dst", dst);
    imshow("src_warp", src_warp);
    
    waitKey(0);
}

結果如下：

視角1的影像視角2的影像視角1的影像校正為視角2

4.2 影像拼接

用 Blender 軟體，獲取相機只做旋轉變換時的視圖1和視圖2，在已知相機標定矩陣和旋轉矩陣的情況下，可計算出兩個視圖之間的單應性矩陣，從而完成拼接，

#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/highgui.hpp"
#include "opencv2/calib3d.hpp"

using namespace cv;

int main()
{
    // 1) read image
    Mat img1 = imread("view1.jpg");
    Mat img2 = imread("view2.jpg");
    if (img1.empty() || img2.empty())
        return -1;

    // 2) camera pose from Blender at location 1
    Mat c1Mo = (Mat_<double>(4, 4) << 0.9659258723258972, 0.2588190734386444, 0.0, 1.5529145002365112,
                                      0.08852133899927139, -0.3303661346435547, -0.9396926164627075, -0.10281121730804443,
                                      -0.24321036040782928, 0.9076734185218811, -0.342020183801651, 6.130080699920654,
                                      0, 0, 0, 1);

    // camera pose from Blender at location 2
    Mat c2Mo = (Mat_<double>(4, 4) << 0.9659258723258972, -0.2588190734386444, 0.0, -1.5529145002365112,
                                      -0.08852133899927139, -0.3303661346435547, -0.9396926164627075, -0.10281121730804443,
                                      0.24321036040782928, 0.9076734185218811, -0.342020183801651, 6.130080699920654,
                                      0, 0, 0, 1);

    // 3) camera intrinsic parameters
    Mat cameraMatrix = (Mat_<double>(3, 3) << 700.0, 0.0, 320.0,
                                              0.0, 700.0, 240.0,
                                              0, 0, 1);
    // 4) extract rotation
    Mat R1 = c1Mo(Range(0, 3), Range(0, 3));
    Mat R2 = c2Mo(Range(0, 3), Range(0, 3));

    // 5) compute rotation displacement: c1Mo * oMc2
    Mat R_2to1 = R1 * R2.t();

    // 6) homography
    Mat H = cameraMatrix * R_2to1 * cameraMatrix.inv();
    H /= H.at<double>(2, 2);

    // 7) warp 
    Mat img_stitch;
    warpPerspective(img2, img_stitch, H, Size(img2.cols * 2, img2.rows));

    // 8) stitch
    Mat half = img_stitch(Rect(0, 0, img1.cols, img1.rows));
    img1.copyTo(half);
    imshow("Panorama stitching", img_stitch);

    waitKey(0);
}

輸出結果如下：

相機視圖1 相機視圖1 拼接后的視圖

4.3 相機位姿估計

一組從不同視角拍攝的標定板，預先知道其拍攝相機的內參，則呼叫 solvePnPRansac() 函式，可估計出該相機拍攝時的位姿

#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/highgui.hpp"
#include "opencv2/calib3d.hpp"

using namespace std;
using namespace cv;

Size kPatternSize = Size(9, 6);
float kSquareSize = 0.025;

int main()
{
    // 1) read image
    Mat src = https://www.cnblogs.com/xinxue/p/imread("images/left01.jpg");
    if (src.empty())
        return -1;
    // prepare for subpixel corner
    Mat src_gray;
    cvtColor(src, src_gray, COLOR_BGR2GRAY);

    // 2) find chessboard corners
    vector<Point2f> corners;
    bool patternfound = findChessboardCorners(src, kPatternSize, corners);

    // 3) get subpixel accuracy
    if (patternfound) {
        cornerSubPix(src_gray, corners, Size(11, 11), Size(-1, -1), TermCriteria(TermCriteria::EPS + TermCriteria::MAX_ITER, 30, 0.1));
    } else {
        return -1;
    }
    
    // 4) define object coordinates
    vector<Point3f> objectPoints;
    for (int i = 0; i < kPatternSize.height; i++) 
    { 
        for (int j = 0; j < kPatternSize.width; j++)
        {
            objectPoints.push_back(Point3f(float(j*kSquareSize), float(i*kSquareSize), 0));
        }
    }

    // 5) input camera intrinsic parameters
    Mat cameraMatrix = (Mat_<double>(3, 3) << 5.3591573396163199e+02, 0.0, 3.4228315473308373e+02, 
                                              0.0, 5.3591573396163199e+02, 2.3557082909788173e+02, 
                                              0.0, 0.0, 1.);
    Mat distCoeffs = (Mat_<double>(5, 1) << -2.6637260909660682e-01, -3.8588898922304653e-02, 1.7831947042852964e-03, 
                                            -2.8122100441115472e-04, 2.3839153080878486e-01);

    // 6) compute rotation and translation vectors
    Mat rvec, tvec;
    solvePnPRansac(objectPoints, corners, cameraMatrix, distCoeffs, rvec, tvec);

    // 7) project estimated pose on the image 
    drawFrameAxes(src, cameraMatrix, distCoeffs, rvec, tvec, 2*kSquareSize);
    imshow("Pose from coplanar points", src);

    waitKey(0);
}