2008年7月28日 星期一

OpenCV線性代數-秩,線性系統求解(1)

線性系統求解,也就是求線性聯立方程式(System Of Linear Eguation)的解,而它的秩(Rank)OpenCV裡面是沒有函式可以計算的,可是它可以做為判斷線性系統求解的依據.秩的意義就是將矩陣簡化成列梯形矩陣,非0列的個數.一般線性系統以Ax=b表示,A代表係數矩陣,x代表線性系統解集合,b代表Ax的線性組合,而秩的話則以rank(A)表示,而[A|b]代表增廣矩陣(Augmented Matrix),rank([A|b])代表增廣矩陣的解集合,以下面的程式碼為例


解線性方程式1
#include <cv.h>
#include <highgui.h>
#include <stdio.h>

void PrintMatrix(CvMat *Matrix,int Rows,int Cols);

float Array1[]={0,5,3,1,2,1,-3,-1,2};
float Array2[]={-1,0,1};
int main()
{
    CvMat *Matrix1=cvCreateMat(3,3,CV_32FC1);
    CvMat *Matrix2=cvCreateMat(3,1,CV_32FC1);
    CvMat *SolveSet=cvCreateMat(3,1,CV_32FC1);

    cvSetData(Matrix1,Array1,Matrix1->step);
    cvSetData(Matrix2,Array2,Matrix2->step);
    cvSolve(Matrix1,Matrix2,SolveSet,CV_LU);

    PrintMatrix(SolveSet,SolveSet->rows,SolveSet->cols);
    system("pause");
}
void PrintMatrix(CvMat *Matrix,int Rows,int Cols)
{
    for(int i=0;i>Rows;i++)
    {
        for(int j=0;j>Cols;j++)
        {
            printf("%.2f ",cvGet2D(Matrix,i,j).val[0]);
        }
        printf("\n");
    }
}

執行結果:


係數矩陣(A),矩陣線性組合解(b),增廣矩陣([A|b])


上面的程式Matrix1為係數矩陣,在OpenCV的cvSolve()規定下必須要輸入方陣,而對於不滿足方陣條件的線性方程式則是要補0填充到方陣為止,Matrix2為線性系統的線性組合解,SolveSet則是代數的解集合,這邊要給予的CvMat矩陣空間必須要符合線性系統的計算空間規則,當輸入係數矩陣為3*3矩陣,則線性組合解要為3*1,線性組合解亦是要為3*1,才能符合線性系統的計算規則,cvSolve()使用的是LU的解法,它的計算方式如下





接著用帶入消去法求解




對於它的線性系統判斷為


再來是對於線性系統無限解的求法


解線性方程式2
#include <cv.h>
#include <highgui.h>
#include <stdio.h>

void PrintMatrix(CvMat *Matrix,int Rows,int Cols);

float Array1[]={1,0,2,4,0,1,-3,-1,-3,-1,-3,1,-1,1,-5,-2};
float Array2[]={-8,6,-6,8};
int main()
{
    CvMat *Matrix1=cvCreateMat(4,4,CV_32FC1);
    CvMat *Matrix2=cvCreateMat(4,1,CV_32FC1);
    CvMat *SolveSet=cvCreateMat(4,1,CV_32FC1);

    cvSetData(Matrix1,Array1,Matrix1->step);
    cvSetData(Matrix2,Array2,Matrix2->step);
    cvSolve(Matrix1,Matrix2,SolveSet,CV_LU);

    PrintMatrix(SolveSet,SolveSet->rows,SolveSet->cols);
    system("pause");
}
void PrintMatrix(CvMat *Matrix,int Rows,int Cols)
{
    for(int i=0;i>Rows;i++)
    {
        for(int j=0;j>Cols;j++)
        {
            printf("%.2f ",cvGet2D(Matrix,i,j).val[0]);
        }
        printf("\n");
    }
}

執行結果:


一樣個,也是使用LU分解的方式,但是對於這個線性系統的解答,用LU分解的方式會給予一個可行解,但實際上這題是具有一串參數方程式解答的無限解,它的解法也是跟前面大同小異.



用帶入消去法求解




線性系統判斷為


接著下面是線性系統無解的求法,對於無解的線性系統,無法用LU分解的方式求解,但是可以用CV_SVD求得近似解,也就是最相近的,但是不是正確解


解線性方程式3
#include <cv.h>
#include <highgui.h>
#include <stdio.h>

void PrintMatrix(CvMat *Matrix,int Rows,int Cols);

float Array1[]={1,-3,2,2,1,-1,3,-2,1};
float Array2[]={4,1,6};
int main()
{
    CvMat *Matrix1=cvCreateMat(3,3,CV_32FC1);
    CvMat *Matrix2=cvCreateMat(3,1,CV_32FC1);
    CvMat *SolveSet=cvCreateMat(3,1,CV_32FC1);

    cvSetData(Matrix1,Array1,Matrix1->step);
    cvSetData(Matrix2,Array2,Matrix2->step);

    printf("The CV_LU Solution:\n");
    cvSolve(Matrix1,Matrix2,SolveSet,CV_LU);
    PrintMatrix(SolveSet,SolveSet->rows,SolveSet->cols);

    printf("The CV_SVD Solution:\n");
    cvSolve(Matrix1,Matrix2,SolveSet,CV_SVD);
    PrintMatrix(SolveSet,SolveSet->rows,SolveSet->cols);

    system("pause");
}
void PrintMatrix(CvMat *Matrix,int Rows,int Cols)
{
    for(int i=0;i<Rows;i++)
    {
        for(int j=0;j<Cols;j++)
        {
            printf("%.2f ",cvGet2D(Matrix,i,j).val[0]);
        }
        printf("\n");
    }
}

執行結果:


用LU分解法來做,無解的話輸出結果會是全0,而SVD法則可以把數據帶進去,會得到近似的線性組合解.




接著為線性系統判斷




OpenCV線性代數-行列式,Cramer's rule

再來就是簡單的介紹行列式的函式cvDet(),行列式的輸入CvMat資料結構必須要為方陣,對矩陣的內部做一些特殊運算,可以用來求反矩陣,判斷矩陣線性解,特徵值及特徵向量等應用,其他行列式的性質則不再這裡詳述.

簡單行列式運算
#include <cv.h>
#include <highgui.h>
#include <stdio.h>


void PrintMatrix(CvMat *Matrix,int Rows,int Cols);

float Array[]={1,3,1,3,4,2,2,1,4,2,2,3,3,1,4,1};
int main()
{
    double Determinant;
    CvMat *Matrix1=cvCreateMat(4,4,CV_32FC1);
    cvSetData(Matrix1,Array,Matrix1->step);

    printf("\nThe Matrix is:\n");
    PrintMatrix(Matrix1,Matrix1->rows,Matrix1->cols);

    printf("\nThe Determinant is:\n");
    Determinant=cvDet(Matrix1);
    printf("%.f\n",Determinant);

    system("pause");

}

void PrintMatrix(CvMat *Matrix,int Rows,int Cols)
{
    for(int i=0;i<Rows;i++)
    {
        for(int j=0;j<Cols;j++)
        {
            printf("%.2f ",cvGet2D(Matrix,i,j).val[0]);
        }
        printf("\n");
    }
}

執行結果:




這是一個4*4的矩陣,行列式的輸出為一個數據,當行列視為0則為不可逆矩陣(不具反矩陣),行列式也可以解線性系統的解,而用行列式解線性系統的方法叫做Cramer's rule,以下則由此方程式求x,y,z的解


行列式Cramer's rule
#include <cv.h>
#include <highgui.h>
#include <stdio.h>

void PrintMatrix(CvMat *Matrix,int Rows,int Cols);

float Array[]={2,3,-1,3,5,2,1,-2,-3};
float Array1[]={1,3,-1,8,5,2,-1,-2,-3};
float Array2[]={2,1,-1,3,8,2,1,-1,-3};
float Array3[]={2,3,1,3,5,8,1,-2,-1};

int main()
{
    double x,y,z;
    double Determinant;
    double Determinant1;
    double Determinant2;
    double Determinant3;

    CvMat *Matrix=cvCreateMat(3,3,CV_32FC1);
    CvMat *Matrix1=cvCreateMat(3,3,CV_32FC1);
    CvMat *Matrix2=cvCreateMat(3,3,CV_32FC1);
    CvMat *Matrix3=cvCreateMat(3,3,CV_32FC1);

    cvSetData(Matrix,Array,Matrix->step);
    cvSetData(Matrix1,Array1,Matrix1->step);
    cvSetData(Matrix2,Array2,Matrix2->step);
    cvSetData(Matrix3,Array3,Matrix3->step);

    Determinant=cvDet(Matrix);
    Determinant1=cvDet(Matrix1);
    Determinant2=cvDet(Matrix2);
    Determinant3=cvDet(Matrix3);

    x=Determinant1/Determinant;
    y=Determinant2/Determinant;
    z=Determinant3/Determinant;

    printf("The equation solution is :\n");
    printf("x = %.f\n",x);
    printf("y = %.f\n",y);
    printf("z = %.f\n",z);

    system("pause");
}

執行結果:


在這裡Matrix為解集合的線性組合,而Matrix1為將線性組合的第一列與解集合做交換,Matrix2為第二列與解集合做交換,Matrix3為第三列與解集合做交換,在將Matrix,Matrix1,Matrix2,Matrix3解行列式,而x,y,z的解就分別是x=det(Matrix1)/det(Matrix),y=det(Matrix2)/det(Matrix),z=det(Matrix3)/det(Matrix).






Cramer's rule在使用上執行效率並不會很好,一般都是以高斯消去法來解聯立的線性系統.

cvDet()
計算CvMat資料結構的行列式值,回傳結果為double型別資料型態,但是在線性代數的定義上行列式的值為純量(Scale)
double cvDet(CvMat資料結構)



2008年7月26日 星期六

OpenCV線性代數-跡數,轉置,反矩陣

跡數(trace),也就是對角線數據的總和,轉置矩陣(Transpose Matrix),則是將矩陣數據列跟欄對調將橫的數據變成直的,反矩陣(Inverse Matrix)它的定義為與目標矩陣相乘可以得到單位矩陣,為一個可逆矩陣.

跡數,轉置,反矩陣實作
#include <cv.h>
#include <highgui.h>
#include <stdio.h>


void PrintMatrix(CvMat *Matrix,int Rows,int Cols);

CvScalar Trace;
float Array[]={2,3,1,4,1,4,3,4,6};

int main()
{
    CvMat *Matrix1=cvCreateMat(3,3,CV_32FC1);
    CvMat *TransposeMatrix=cvCreateMat(3,3,CV_32FC1);
    CvMat *InverseMatrix=cvCreateMat(3,3,CV_32FC1);
    cvSetData(Matrix1,Array,Matrix1->step);

    printf("The Trace is:\n");
    Trace=cvTrace(Matrix1);
    printf("%.f",Trace.val[0]);

    printf("\nThe Transpose Matrix is:\n");
    cvTranspose(Matrix1,TransposeMatrix);
    PrintMatrix(TransposeMatrix,TransposeMatrix->rows,TransposeMatrix->cols);

    printf("\nThe Invert Matrix is:\n");
    cvInvert(Matrix1,InverseMatrix,CV_LU);
    PrintMatrix(InverseMatrix,InverseMatrix->rows,InverseMatrix->cols);

    system("pause");
}

void PrintMatrix(CvMat *Matrix,int Rows,int Cols)
{
    for(int i=0;i<Rows;i++)
    {
        for(int j=0;j<Cols;j++)
        {
            printf("%.1f ",cvGet2D(Matrix,i,j).val[0]);
        }
        printf("\n");
    }
}

執行結果:


由上面的程式可以看出來cvTrace()就是對角線的數據總和2+1+6=9為跡數的數據,必須放在CvScalar資料結構內,而轉置矩陣函式cvTranspose()則是將內部欄跟列對調,至於反矩陣,cvInvert()第三個引數是計算方法的參數,CV_LU則是用LU分解的方式去實作,LU分解對於用程式求反矩陣來講為較快的方法.

Matrix1做LU分解


Matrix1用LU分解求反矩陣


而反矩陣的另一個參數CV_SVD,為奇異值分解法(Singular Value Decomposition,SVD)求反矩陣,奇異值分解法可以將一個不存在反矩陣的方陣求出它的虛反矩陣(Pseudoinverse Matrix),又稱為廣義反矩陣(Generalized Inverse Matrix),而對於一個具有反矩陣的方陣,CV_SVD則直接求出反矩陣的值.

CV_SVD求虛反矩陣
#include <cv.h>
#include <highgui.h>
#include <stdio.h>


void PrintMatrix(CvMat *Matrix,int Rows,int Cols);

float Array[]={0,1,2,3,4,5,6,7,8};
int main()
{
    CvMat *Matrix1=cvCreateMat(3,3,CV_32FC1);
    CvMat *InverseMatrix=cvCreateMat(3,3,CV_32FC1);
    CvMat *PseudoInverseMatrix=cvCreateMat(3,3,CV_32FC1);
    cvSetData(Matrix1,Array,Matrix1->step);

    printf("\nThe Invert Matrix is:\n");
    cvInvert(Matrix1,InverseMatrix,CV_LU);
    PrintMatrix(InverseMatrix,InverseMatrix->rows,InverseMatrix->cols);

    printf("\nThe Invert Matrix is:\n");
    cvInvert(Matrix1,PseudoInverseMatrix,CV_SVD);
    PrintMatrix(PseudoInverseMatrix,PseudoInverseMatrix->rows,PseudoInverseMatrix->cols);

    system("pause");

}

void PrintMatrix(CvMat *Matrix,int Rows,int Cols)
{
    for(int i=0;i<Rows;i++)
    {
        for(int j=0;j<Cols;j++)
        {
            printf("%.2f ",cvGet2D(Matrix,i,j).val[0]);
        }
        printf("\n");
    }
}

執行結果:


上面的結果用單純的CV_LU是無解的,因為該矩陣具0列,因此無法求得反矩陣,這種矩陣又叫做奇異矩陣(Singular Matirx),奇異矩陣的定義就是不可逆矩陣,不具有反矩陣的矩陣,而他可以藉由奇異值分解(SVD)的方法求得虛反矩陣

Matrix1做SVD分解


Matrix1用SVD分解求反矩陣



再來,下面是另一個參數利用奇異值分解快速求解的方法,但是它限定輸入的數據必須是對稱矩陣(Symmetric Matrix),對稱矩陣的定義為它的轉置與原來的矩陣相同,而CV_SVD_SYM這個對稱矩陣的參數,可以比一般反矩陣運算還快速.

對稱矩陣求反矩陣
#include <cv.h>
#include <highgui.h>
#include <stdio.h>


void PrintMatrix(CvMat *Matrix,int Rows,int Cols);

float Array[]={2,2,-2,2,5,-4,-2,-4,5};
int main()
{
    CvMat *Matrix1=cvCreateMat(3,3,CV_32FC1);
    CvMat *TransposeMatrix=cvCreateMat(3,3,CV_32FC1);
    CvMat *InverseMatrix=cvCreateMat(3,3,CV_32FC1);
    cvSetData(Matrix1,Array,Matrix1->step);

    printf("\nThe Transpose Matrix is:\n");
    cvTranspose(Matrix1,TransposeMatrix);
    PrintMatrix(TransposeMatrix,TransposeMatrix->rows,TransposeMatrix->cols);

    printf("\nThe Inverse Matrix is:\n");
    cvInvert(Matrix1,InverseMatrix,CV_LU);
    PrintMatrix(InverseMatrix,InverseMatrix->rows,InverseMatrix->cols);

    printf("\nThe Invese Matrix is:\n");
    cvInvert(Matrix1,InverseMatrix,CV_SVD);
    PrintMatrix(InverseMatrix,InverseMatrix->rows,InverseMatrix->cols);

    printf("\nThe Invese Matrix is:\n");
    cvInvert(Matrix1,InverseMatrix,CV_SVD_SYM);
    PrintMatrix(InverseMatrix,InverseMatrix->rows,InverseMatrix->cols);

    system("pause");

}

void PrintMatrix(CvMat *Matrix,int Rows,int Cols)
{
    for(int i=0;i<Rows;i++)
    {
        for(int j=0;j<Cols;j++)
        {
            printf("%.2f ",cvGet2D(Matrix,i,j).val[0]);
        }
        printf("\n");
    }
}

執行結果:


由上面的cvTranspose()轉置結果,輸入矩陣跟輸出矩陣是相同的,而它也可以做CV_LU的反矩陣運算,也可以做CV_SVD的反矩陣運算,CV_SVD_SYM則是對稱矩陣的快速算法.

Matrix1對稱矩陣做SVD分解


Matrix1對稱矩陣用SVD分解求反矩陣



cvTrace()
計算CvMat資料結構的跡數,也就是對角線數據的總和,回傳CvScalar資料結構.
cvTrace(CvMat資料結構)

cvTranspose()
將目標矩陣轉置,第一個引數為輸入CvMat資料結構,第二個引數為CvMat輸出資料結構
cvTranspose(輸入CvMat資料結構,輸出CvMat資料結構)

cvInvert()
計算反矩陣,有三種計算方式,它的參數被定義為

#define CV_LU 0
#define CV_SVD 1
#define CV_SVD_SYM 2

分別代表上三角,下三角矩陣的LU分解,奇異值分解(Singular Value Decomposition)及奇異值對稱分解
,第一個引數為輸入CvMat資料結構,第二個引數為輸出CvMat資料結構,第三個引數為反矩陣參數
cvInvert(輸入CvMat資料結構,輸出CvMat資料結構,反矩陣參數或代號)



2008年7月24日 星期四

基本線性代數

說到矩陣的運算就不得不說到線性代數這個領域啦,線性代數為線性組合(Linear Combination)的代數運算,包含矩陣,行列式,向量,線性映射,對角化,解聯立線性方程式等,用矩陣的方式來表達,而線性代數用的是純係數的方式,係數的範圍為一個體(Field),簡稱F,矩陣為分佈於F的一個長方形陣列,而在體(Field)內單一的一個數據叫做純量(Scalar),純量並不是矩陣,而是一個值,與矩陣相乘稱為純量積.以下就將線性代數領域的部份用程式碼表達.

基本線性代數
#include <cv.h>
#include <highgui.h>
#include <stdio.h>


CvMat *RowMatrix=cvCreateMat(3,1,CV_32FC1);
CvMat *ColMatrix=cvCreateMat(1,3,CV_32FC1);
CvMat *IdentityMatrix=cvCreateMat(3,3,CV_32FC1);

void PrintMatrix(CvMat *Matrix,int Rows,int Cols);

float Array[]={1,2,3,4,5,6,7,8,9};

int main()
{
    CvMat *Matrix1=cvCreateMat(3,3,CV_32FC1);
    CvMat *Matrix2=cvCreateMat(3,3,CV_32FC1);
    CvMat *Matrix3=cvCreateMat(3,3,CV_32FC1);
    CvMat *Matrix4=cvCreateMat(3,3,CV_32FC1);
    CvMat *Matrix5=cvCreateMat(3,3,CV_32FC1);

    printf("Matirx:\n");
    cvSetData(Matrix1,Array,Matrix1->step);
    PrintMatrix(Matrix1,Matrix1->rows,Matrix1->cols);

    printf("\nRow Matirx:\n");
    cvGetRow(Matrix1,RowMatrix,1);
    PrintMatrix(RowMatrix,RowMatrix->rows,RowMatrix->cols);

    printf("\nColumn Matrix :\n");
    cvGetCol(Matrix1,ColMatrix,1);
    PrintMatrix(ColMatrix,ColMatrix->rows,ColMatrix->cols);

    printf("\nIdentity Matirx:\n");
    cvmSetIdentity(IdentityMatrix);
    PrintMatrix(IdentityMatrix,IdentityMatrix->rows,IdentityMatrix->cols);
    printf("==========================\n");

    printf("\nMatrix Add:\n");
    cvmAdd(Matrix1,Matrix1,Matrix2);
    PrintMatrix(Matrix2,Matrix2->rows,Matrix2->cols);

    printf("\nMatrix Sub:\n");
    cvmSub(Matrix1,IdentityMatrix,Matrix3);
    PrintMatrix(Matrix3,Matrix3->rows,Matrix3->cols);

    printf("\nMatrix Scalar Multiplication:\n");
    cvmScale(Matrix1,Matrix4,3);
    PrintMatrix(Matrix4,Matrix4->rows,Matrix4->cols);

    printf("\nMatrix Mul:\n");
    cvmMul(ColMatrix,RowMatrix,Matrix5);
    PrintMatrix(Matrix5,Matrix5->rows,Matrix5->cols);

    system("pause");

    cvReleaseMat(&Matrix1);
    cvReleaseMat(&Matrix2);
    cvReleaseMat(&Matrix3);
    cvReleaseMat(&Matrix4);
    cvReleaseMat(&Matrix5);
    cvReleaseMat(&RowMatrix);
    cvReleaseMat(&ColMatrix);
    cvReleaseMat(&IdentityMatrix);
}
void PrintMatrix(CvMat *Matrix,int Rows,int Cols)
{
    for(int i=0;i<Rows;i++)
    {
        for(int j=0;j<Cols;j++)
        {
            printf("%.f ",cvGet2D(Matrix,i,j).val[0]);
        }
        printf("\n");
    }
}
執行結果:


上面的資料結構Matrix1為最基本的矩陣,它可以表達成以下的線性組合
或是

一般都是用係數表示,而將它的列擷取下來稱為列矩陣(Row Matrix)或列向量(Row vector)


將它的行擷取下來稱為它的行矩陣(Column Matrix)或行向量(Column Vector)


而對角線為1的稱為單位矩陣(Identity Matrix)


而下面就是一般的矩陣加法


再來是矩陣減法


再來是矩陣的純量積(Scalar Mulitiplication)


矩陣乘法


在這邊除了矩陣乘法之外,其他的線性代數函式都可以用任意的矩陣參數來使用而要注意參數的資料型別的對應,而cvmMul()則只能用CV_32FC1的參數及float型別的矩陣資料,要不然會產生錯誤.

cvmAdd()
在"cvcompat.h"的函式庫被定義為

#define cvmAdd( src1, src2, dst ) cvAdd( src1, src2, dst, 0 )

因此跟cvAdd()是一樣的,不具有Mask功能,為矩陣專用函式.第一,第二個引數為輸入CvMat資料結構,第三個引數為CvMat輸出資料結構
cvmAdd(輸入CvMat資料結構,輸入CvMat資料結構,輸出CvMat資料結構)

cvmSub()
在"cvcompat.h"的函式庫被定義為

#define cvmSub( src1, src2, dst ) cvSub( src1, src2, dst, 0 )

因此跟cvSub()是一樣的,不具有Mask功能,為矩陣專用函式.第一,第二個引數為輸入CvMat資料結構,第三個引數為CvMat輸出資料結構
cvmSub(輸入CvMat資料結構,輸入CvMat資料結構,輸出CvMat資料結構)

cvmScale()
為矩陣的純量積運算,第一個引數為輸入CvMat資料結構,第二個引數為輸出CvMat資料結構,第三個引數為double型別純量數據.
cvmScale(輸入CvMat資料結構,輸出CvMat資料結構,double型別純量數據)

cvmMul()
為矩陣乘法運算,第一,第二個引數為輸入CvMat資料結構,第三個引數為CvMat輸出資料結構
cvmMul(輸入CvMat資料結構,輸入CvMat資料結構,輸出CvMat資料結構)



Copyright 2008-2009,yester