amatrix2.cpp

来自「微软的基于HMM的人脸识别原代码, 非常经典的说」· C++ 代码 · 共 540 行 · 第 1/2 页

/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                        Intel License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2000, Intel Corporation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of Intel Corporation may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#define DOT     0
#define CROSS   1
#define TRACE   2
#define MT      3
#define PP      4
#define EV      5
#define DET     6

#include "CvmTest.h"
#include <float.h>
#define PI  3.141593f
#define FEPS 1.0e-5
#define DEPS 1.0e-10

/* Testing parameters */
static char* func_names[] =
{
    "cvmDotProduct",
    "cvmCrossProduct",
    "cvmTrace",
    "cvmMulTransposed",
    "cvmPerspectiveProject",
	"cvmEigenVV",
	"cvmDet"
};

static char* test_desc = "Matrix - Matrix operations (2)";
static int max_length, length1, length2, ntests;
static int step;
static int read_param = 0;
static int test_type=6;

#include "_cvMatrixTestFunctions.cpp"

static int fcaMatrMatr2( void* prm )
{
    int i, j, it;
    const int nh = 20;
    int        Errors = 0;
    long lParam = (long)prm;
    double result, result0;
    CvMat fSrc1, fSrc2, fDst1, fDst2, dSrc1, dSrc2, dDst1, dDst2, dDst1a, dDst2a;
    CvMatArray fSrc, dSrc, fDst, dDst;
    AtsRandState state, stt;
      
    /* Initialization global parameters */
    if( !read_param )
    {
        /* Determine which tests are needed to run */
        trsCaseRead( &test_type, "/D/C/T/M/P/E/DET/a", "a",
        "a - all, D - dot product, C - cross product, T - trace, M - multiply transpose, "
        "P - perspective projection, E - eigen values and vectors, DET - determinant" );
        
        /* Reading test-parameters */
        trsiRead( &ntests,     "25", " number of tests" );
        trsiRead( &max_length, "32", " maximal size of matrices" );
        read_param = 1;
    }
    if((test_type!=7)&&(test_type!=lParam)) return TRS_UNDEF;

    atsRandInit( &state, -256.0, 256.0, 10 );
    atsRandInit( &stt, 1.001, max_length+0.999, 13 );

    for( it=0; it<ntests; it++ )
    {

        switch( lParam )
        {
		case DOT:
                length1 = atsRand32s( &stt );

                fSrc1 = cvMat(length1, 1, CV_MAT32F, NULL);
                fSrc2 = cvMat(length1, 1, CV_MAT32F, NULL);
                dSrc1 = cvMat(length1, 1, CV_MAT64D, NULL);
                dSrc2 = cvMat(length1, 1, CV_MAT64D, NULL);
            	cvmAlloc(&fSrc1);
            	cvmAlloc(&fSrc2);
            	cvmAlloc(&dSrc1);
            	cvmAlloc(&dSrc2);

                atsbRand32f( &state, fSrc1.data.fl, length1 );
                atsbRand32f( &state, fSrc2.data.fl, length1 );
                __fl2db( fSrc1.data.fl, dSrc1.data.db, length1 );
                __fl2db( fSrc2.data.fl, dSrc2.data.db, length1 );
                result  = cvmDotProduct( &fSrc1, &fSrc2 );
                result0 =  __DotProduct( dSrc1.data.db, dSrc2.data.db, length1 );
                if( fabs( ( result - result0 )/result0 ) > FEPS ) Errors++;
        
                atsbRand64d( &state, dSrc1.data.db, length1 );
                atsbRand64d( &state, dSrc2.data.db, length1 );
                result  = cvmDotProduct( &dSrc1, &dSrc2 );
                result0 =  __DotProduct( dSrc1.data.db, dSrc2.data.db, length1 );
                if( fabs( ( result - result0 )/result0 ) > DEPS ) Errors++;

                cvmFree(&fSrc1);
                cvmFree(&fSrc2);
                cvmFree(&dSrc1);
                cvmFree(&dSrc2);
            break;

		case CROSS:
                fSrc1 = cvMat(3, 1, CV_MAT3x1_32F, NULL);
                fSrc2 = cvMat(3, 1, CV_MAT3x1_32F, NULL);
                dSrc1 = cvMat(3, 1, CV_MAT3x1_64D, NULL);
                dSrc2 = cvMat(3, 1, CV_MAT3x1_64D, NULL);
                fDst1 = cvMat(3, 1, CV_MAT3x1_32F, NULL);
                fDst2 = cvMat(3, 1, CV_MAT3x1_32F, NULL);
                dDst1 = cvMat(3, 1, CV_MAT3x1_64D, NULL);
                dDst2 = cvMat(3, 1, CV_MAT3x1_64D, NULL);
            	cvmAlloc(&fSrc1);
            	cvmAlloc(&fSrc2);
            	cvmAlloc(&dSrc1);
            	cvmAlloc(&dSrc2);
            	cvmAlloc(&fDst1);
            	cvmAlloc(&fDst2);
            	cvmAlloc(&dDst1);
            	cvmAlloc(&dDst2);

                atsbRand32f( &state, fSrc1.data.fl, 3 );
                atsbRand32f( &state, fSrc2.data.fl, 3 );
                __fl2db( fSrc1.data.fl, dSrc1.data.db, 3 );
                __fl2db( fSrc2.data.fl, dSrc2.data.db, 3 );
                __CrossProduct( dSrc1.data.db, dSrc2.data.db, dDst2.data.db );
                cvmCrossProduct( &fSrc1, &fSrc2, &fDst1 );
                __fl2db( fDst1.data.fl, dDst1.data.db, 3 );
                Errors += atsCompDoublePrec( dDst1.data.db ,dDst2.data.db, 3, FEPS );
        
                cvmFree(&fSrc1);
                cvmFree(&fSrc2);
                cvmFree(&dSrc1);
                cvmFree(&dSrc2);
                cvmFree(&fDst1);
                cvmFree(&fDst2);
                cvmFree(&dDst1);
                cvmFree(&dDst2);
            break;

        case TRACE:
                length1 = atsRand32s( &stt );
                length2 = atsRand32s( &stt );

                fSrc1 = cvMat(length1, length2, CV_MAT32F, NULL);
                dSrc1 = cvMat(length1, length2, CV_MAT64D, NULL);
            	cvmAlloc(&fSrc1);
            	cvmAlloc(&dSrc1);

                atsbRand32f( &state, fSrc1.data.fl, length1*length2 );
                __fl2db( fSrc1.data.fl, dSrc1.data.db, length1*length2 );
                result  = cvmTrace( &fSrc1 );
                result0 =  __Trace( dSrc1.data.db, length1, length2 );
                if( fabs( ( result - result0 )/result0 ) > FEPS ) Errors++;
        
                atsbRand64d( &state, dSrc1.data.db, length1*length2 );
                result  = cvmTrace( &dSrc1 );
                result0 =  __Trace( dSrc1.data.db, length1, length2 );
                if( fabs( ( result - result0 )/result0 ) > DEPS ) Errors++;

                cvmFree(&fSrc1);
                cvmFree(&dSrc1);
            break;

        case MT:
                length1 = atsRand32s( &stt );
                length2 = atsRand32s( &stt );

                fSrc1 = cvMat(length1, length2, CV_MAT32F, NULL);
                fDst1 = cvMat(length1, length1, CV_MAT32F, NULL);
                fDst2 = cvMat(length2, length2, CV_MAT32F, NULL);
                dSrc1 = cvMat(length1, length2, CV_MAT64D, NULL);
                dDst1 = cvMat(length1, length1, CV_MAT64D, NULL);
                dDst2 = cvMat(length2, length2, CV_MAT64D, NULL);
                dDst1a= cvMat(length1, length1, CV_MAT64D, NULL);
                dDst2a= cvMat(length2, length2, CV_MAT64D, NULL);
            	cvmAlloc(&fSrc1);
            	cvmAlloc(&dSrc1);
            	cvmAlloc(&fDst1);
            	cvmAlloc(&fDst2);
            	cvmAlloc(&dDst1);
            	cvmAlloc(&dDst2);
            	cvmAlloc(&dDst1a);
            	cvmAlloc(&dDst2a);

                atsbRand32f( &state, fSrc1.data.fl, length1*length2 );
                __fl2db( fSrc1.data.fl, dSrc1.data.db, length1*length2 );
                cvmMulTransposed( &fSrc1, &fDst1, 0 );
                __MulTrans( dSrc1.data.db, dDst1.data.db, length1, length2, 0 );
                for( i=0; i<length1*length1; i++ )
                {
                    double r1 = (double)fDst1.data.fl[i];
                    double r2 =         dDst1.data.db[i];
                    if( fabs( ( r1 - r2 )/r2 ) > 50.0*FEPS ) { Errors++;/* printf( " fl., var. 1:  %g  %g  err = %g\n", r1, r2, fabs( ( r1 - r2 )/r2 ) );*/ }
                }
        
                cvmMulTransposed( &fSrc1, &fDst2, 1 );
                __MulTrans( dSrc1.data.db, dDst2.data.db, length1, length2, 1 );
                for( i=0; i<length2*length2; i++ )
                {
                    double r1 = (double)fDst2.data.fl[i];
                    double r2 =         dDst2.data.db[i];
                    if( fabs( ( r1 - r2 )/r2 ) > 50.0*FEPS ) { Errors++;/* printf( " fl., var. 2:  %g  %g  err = %g\n", r1, r2, fabs( ( r1 - r2 )/r2 ) );*/ }
                }
        
                atsbRand64d( &state, dSrc1.data.db, length1*length2 );
                cvmMulTransposed( &dSrc1, &dDst1a, 0 );
                __MulTrans( dSrc1.data.db, dDst1.data.db, length1, length2, 0 );
                for( i=0; i<length1*length1; i++ )
                {
                    double r1 = dDst1a.data.db[i];
                    double r2 =  dDst1.data.db[i];
                    if( fabs( ( r1 - r2 )/r2 ) > DEPS ) { Errors++;/* printf( " db., var. 1:  %g  %g  err = %g\n", r1, r2, fabs( ( r1 - r2 )/r2 ) ); */}
                }
        
                cvmMulTransposed( &dSrc1, &dDst2a, 1 );
                __MulTrans( dSrc1.data.db, dDst2.data.db, length1, length2, 1 );
                for( i=0; i<length2*length2; i++ )
                {
                    double r1 = dDst2a.data.db[i];
                    double r2 =  dDst2.data.db[i];
                    if( fabs( ( r1 - r2 )/r2 ) > DEPS ) { Errors++;/* printf( " db., var. 2:  %g  %g  err = %g\n", r1, r2, fabs( ( r1 - r2 )/r2 ) ); */}
                }

                cvmFree(&fSrc1);
                cvmFree(&dSrc1);
                cvmFree(&fDst1);
                cvmFree(&dDst1);
                cvmFree(&fDst2);
                cvmFree(&dDst2);
                cvmFree(&dDst1a);
                cvmFree(&dDst2a);
            break;