perspective transform透視矩陣快速求法+矩形矯正

算了半天一直在思考如何快速把矩陣算出來，網上基本都是在說邊長為1的正方形的變換方式=。=   不怎么用得上…… 公式推導推半天，計算還麻煩。。。。

 ++++++++++++++++++++++++++++++

對於透視變換：

其中我們在使用的時候最重要的是算出系數……………………  這里我們按照i=1計算

找了半天找到了這一篇論文：http://www.ixueshu.com/document/027c165d22e18077318947a18e7f9386.html  

 

 

里面提到我們可以將公式

 

改寫成的模式。

那么根據矩形矯正前后的8個頂點可以得到公式：

那么只要計算中間8x8矩陣的逆，乘上變換后的uv矩陣，就可以求出所有待定系數。其中點的檢測需要用到hough算法，點和點之間的對應我做了個交互，不然太麻煩了。。。。

 

矩陣求逆代碼：其中N為8

#include <iostream>
using namespace std;
#define N 8    //測試矩陣維數定義

//按第一行展開計算|A|
double getA(double arcs[N][N], int n)
{
    if (n == 1)
    {
        return arcs[0][0];
    }
    double ans = 0;
    double temp[N][N] = { 0.0 };
    int i, j, k;
    for (i = 0; i<n; i++)
    {
        for (j = 0; j<n - 1; j++)
        {
            for (k = 0; k<n - 1; k++)
            {
                temp[j][k] = arcs[j + 1][(k >= i) ? k + 1 : k];

            }
        }
        double t = getA(temp, n - 1);
        if (i % 2 == 0)
        {
            ans += arcs[0][i] * t;
        }
        else
        {
            ans -= arcs[0][i] * t;
        }
    }
    return ans;
}

//計算每一行每一列的每個元素所對應的余子式，組成A*
void  getAStart(double arcs[N][N], int n, double ans[N][N])
{
    if (n == 1)
    {
        ans[0][0] = 1;
        return;
    }
    int i, j, k, t;
    double temp[N][N];
    for (i = 0; i<n; i++)
    {
        for (j = 0; j<n; j++)
        {
            for (k = 0; k<n - 1; k++)
            {
                for (t = 0; t<n - 1; t++)
                {
                    temp[k][t] = arcs[k >= i ? k + 1 : k][t >= j ? t + 1 : t];
                }
            }


            ans[j][i] = getA(temp, n - 1);  //此處順便進行了轉置
            if ((i + j) % 2 == 1)
            {
                ans[j][i] = -ans[j][i];
            }
        }
    }
}

//得到給定矩陣src的逆矩陣保存到des中。
bool GetMatrixInverse(double src[N][N], int n, double des[N][N])
{
    double flag = getA(src, n);
    double t[N][N];
    if (0 == flag)
    {
        cout << "原矩陣行列式為0，無法求逆。請重新運行" << endl;
        return false;//如果算出矩陣的行列式為0，則不往下進行
    }
    else
    {
        getAStart(src, n, t);
        for (int i = 0; i<n; i++)
        {
            for (int j = 0; j<n; j++)
            {
                des[i][j] = t[i][j] / flag;
            }

        }
    }

    return true;
}

 

 主函數內代碼：（僅部分）

vector<pair<int, int>> P ,newPoint; //P 變換后標准A4的點  newPoint 變換前的點   計算時候計算變換后變換到變換前

//省略中間數值插入和輸出圖片構成

P.push_back(make_pair(0, 0));
    P.push_back(make_pair((int)len1, 0));
    P.push_back(make_pair((int)len1, (int)len2));
    P.push_back(make_pair(0, (int)len2));

    int uv[8] = { newPoint[0].first, newPoint[0].second,
        newPoint[1].first, newPoint[1].second,
        newPoint[2].first, newPoint[2].second,
        newPoint[3].first, newPoint[3].second };
    
    double src[8][8] =
    { { P[0].first, P[0].second, 1, 0, 0, 0, -newPoint[0].first*P[0].first, -newPoint[0].first*P[0].second },
    { 0, 0, 0, P[0].first, P[0].second, 1, -newPoint[0].second*P[0].first, -newPoint[0].second*P[0].second },

    { P[1].first, P[1].second, 1, 0, 0, 0, -newPoint[1].first*P[1].first, -newPoint[1].first*P[1].second },
    { 0, 0, 0, P[1].first, P[1].second, 1, -newPoint[1].second*P[1].first, -newPoint[1].second*P[1].second },

    { P[2].first, P[2].second, 1, 0, 0, 0, -newPoint[2].first*P[2].first, -newPoint[2].first*P[2].second },
    { 0, 0, 0, P[2].first, P[2].second, 1, -newPoint[2].second*P[2].first, -newPoint[2].second*P[2].second },

    { P[3].first, P[3].second, 1, 0, 0, 0, -newPoint[3].first*P[3].first, -newPoint[3].first*P[3].second },
    { 0, 0, 0, P[3].first, P[3].second, 1, -newPoint[3].second*P[3].first, -newPoint[3].second*P[3].second } };



    double matrix_after[N][N]{};
    bool flag = GetMatrixInverse(src, N, matrix_after);
    if (false == flag) {
        cout << "求不出系數" << endl;
        return;
    }


    cout << "逆矩陣：" << endl;

    for (int i = 0; i<8; i++)
    {
        for (int j = 0; j<8; j++)
        {
            cout << matrix_after[i][j] << " ";
            //cout << *(*(matrix_after+i)+j)<<" ";
        }
        cout << endl;
    }

    double xs[8];
    for (int i = 0; i < 8; i++) {
        double sum = 0;
        for (int t = 0; t < 8; t++) {
            sum += matrix_after[i][t] * uv[t];
        }
        xs[i] = sum;
    }
        //矩形矯正  沒有用雙線性插值
    cimg_forXY(outputimg, x, y) {
        double px = xs[0] * x + xs[1] * y + xs[2];
        double py = xs[3] * x + xs[4] * y + xs[5];
        double p = xs[6] * x + xs[7] * y + 1;

        int u = px / p;
        int v = py / p;

        outputimg(x, y, 0) = paint(u, v, 0);
        outputimg(x, y, 1) = paint(u, v, 1);
        outputimg(x, y, 2) = paint(u, v, 2);
    }

 

以上~

基本實驗了一下：

結果：

 

還可以吧。。。