2008年8月3日 星期日

OpenCV線性代數-cvEigenVV實作

cvEigenVV()為計算CvMat方陣的特徵值(Eigen Value)跟特徵向量(Eigen Vector)的函式,但是,OpenCV的Eigen Value跟Eigen Vector計算並不是一般性的用法,無法處理正常的特徵值跟特徵向量的計算,cvEigenVV()是使用到Jacobi Eigenvalue Algorithm的方法,輸入必須要對稱矩陣,輸入的對稱矩陣做Jacobi Transformation的轉換,這部份的資料就要參考數值分析(Numerical Recipes)的相關書籍了,而Jacobi Eigenvalue Algorithm有將Eigen Value重新Sort過,因此特徵值的排序會是由大到小排列,而不適用於很多矩陣對角化的解題應用.


cvEigenVV()實作
#include <cv.h>
#include <highgui.h>
#include <stdio.h>

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

double Array1[]={2,3,0,3,5,-1,0,-1,2};

int main()
{
    CvMat *Matrix1=cvCreateMat(3,3,CV_64FC1);
    CvMat *EigenValue_Row=cvCreateMat(3,1,CV_64FC1);
    CvMat *EigenValue=cvCreateMat(3,3,CV_64FC1);
    CvMat *EigenVector=cvCreateMat(3,3,CV_64FC1);
    CvMat *EigenVector_Invert=cvCreateMat(3,3,CV_64FC1);
    CvMat *ResultMatrix=cvCreateMat(3,3,CV_64FC1);
    cvSetData(Matrix1,Array1,Matrix1->step);
    cvSetZero(EigenValue);

    cvEigenVV(Matrix1,EigenVector,EigenValue_Row,DBL_EPSILON);

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

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

    printf("\nThe EigenValue is:\n");
    cvSet2D(EigenValue,0,0,cvGet2D(EigenValue_Row,0,0));
    cvSet2D(EigenValue,1,1,cvGet2D(EigenValue_Row,1,0));
    cvSet2D(EigenValue,2,2,cvGet2D(EigenValue_Row,2,0));
    PrintMatrix(EigenValue,EigenValue->rows,EigenValue->cols);

    cvTranspose(EigenVector,EigenVector);
    cvmMul(EigenVector,EigenValue,ResultMatrix);
    cvInvert(EigenVector,EigenVector_Invert,CV_LU);
    cvmMul(ResultMatrix,EigenVector_Invert,ResultMatrix);

    printf("\nTo validate Matrix\n");
    PrintMatrix(ResultMatrix,ResultMatrix->rows,ResultMatrix->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");
    }
}
執行結果:


由上面可以知道Eigen Value的排列方式為


而Eigen Vector的排列方式為


因此,要驗證是否可行,利用


的方法將它還原為對稱矩陣,則必須要將輸出的列矩陣補0讓它成為對角矩陣,並且將Eigen Vector矩陣做轉置,帶入驗證公式,輸出結果為原對稱矩陣,因此可以證明此函式及輸入值無誤

cvEigenVV()第一的引數為特徵向量(Eigen Vector)第二個引數為特徵值(Eigen Value)第三個引數為精確度,DBL_EPSILON為最小double型別精確度,而一般精確度的使用定義為

#define FLT_EPSILON 1.19209290E-07F
#define DBL_EPSILON 2.2204460492503131E-16

這可以來做浮點數的精確度比較運算使用.
而他的使用方式就如同下所述

bool IsEqual(double x,double y)
{
    if(x-y<DBL_EPSILON)
    {
        return true;
    }
    else
    {
        return false;
    }
}

double型別則對照DBL_EPSILON參數,float型別則對照FLT_EPSILON,而cvEigenVV()則可以自行定義它的精確度比較運算的大小.

cvEigenVV()
利用Jacobi Eigenvalue Algorithm法計算Eigen Value,Eigen Vector,輸入只能為對稱矩陣,輸出則是Eigen Value的列舉陣及Eigen Vector的列舉陣,第一個引數為要計算的CvMat資料結構對稱矩陣,第二個引數為CvMat資料結構的Eigen Value列舉陣,第三個引數是CvMat資料結構的Eigen Vector列舉陣,第四個引數為精確度比較運算的大小.
cvEigenVV(輸入CvMat資料結構對稱矩陣,輸出CvMat資料結構的EigenValue列舉陣,輸出CvMat資料結構的EigenVector列舉陣,浮點型別精準度誤差比較數值)



2008年8月1日 星期五

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

在cvSolve()裡面,輸入的係數矩陣規定要方陣,那假如所求的線性系統不是方陣怎麼辦呢?那當然就是要用補0了方式填充成方陣,使用的方法請看以下程式碼.


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

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

float Array1[]={1,1,0,2,1,0,3,2,0};
float Array2[]={1,3,4};
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_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");
    }
}

執行結果:



線性系統判斷


上面的程式是尾端補0,因為這個線性系統只存在兩個未知數,三個方程式,要給它弄成方陣,則只好假設有第三個未知數存在,但係數為0,而只要是需要補0的線性系統運算,全部都需要用CV_SVD來解,LU是無法解的.再來,下面是整欄補零的範例.


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

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

float Array1[]={1,1,1,1,2,0,0,0,0};
float Array2[]={4,6,0};
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_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");
    }
}

執行結果:


線性系統判斷


這個線性系統是四個未知數,三個方程式,很顯然的是無限解,用SVD的方法的話就會任意給一個可行解,而至於cvSolve()的另外一個參數CV_SVD_SYM,則是線性系統對對稱係數矩陣的加速算法,也是跟前面相關的程式碼同樣的用法,這邊不再累述.

cvSolve()
為解線性方程式的函式,可以對唯一解,無解,無限解,而無限解的問題會給予一個答案,無解的問題則是用CV_SVD參數則可得到近似解,而無限解的問題則用cvSolve()參數可以得到可行解.第一個引數為CvMat係數矩陣,第二個引數為線性方程式的線性組合解,第三個引數為代數的解集合,第四個為解線性方程的參數輸入,分別為CV_LU,CV_SVD,CV_SVD_SYM,與求反矩陣的函式相同參數.
cvSolve(CvMat係數矩陣結構,CvMat結構線性組合解,CvMat結構代數解集合,輸入計算參數或代號)



Copyright 2008-2009,yester