cxmatrix.cpp

opencv库在TI DM6437上的移植,目前包括两个库cv.lib和cxcore.lib的工程

    }
    else
    {
        if( !CV_IS_MAT(srcB))
            CV_CALL( srcB = cvGetMat( srcB, &stubB ));

        if( !CV_IS_MAT(dst))
            CV_CALL( dst = cvGetMat( dst, &dstub ));

        if( !CV_ARE_TYPES_EQ( srcA, srcB ) ||
            !CV_ARE_TYPES_EQ( srcB, dst ))
            CV_ERROR( CV_StsUnmatchedFormats, "" );
    }

    if( !CV_ARE_SIZES_EQ( srcA, srcB ) || !CV_ARE_SIZES_EQ( srcB, dst ))
        CV_ERROR( CV_StsUnmatchedSizes, "" );

    if( CV_MAT_DEPTH(type) == CV_32F )
    {
        float* dstdata = (float*)(dst->data.ptr);
        const float* src1data = (float*)(srcA->data.ptr);
        const float* src2data = (float*)(srcB->data.ptr);

        if( CV_IS_MAT_CONT(srcA->type & srcB->type & dst->type) )
        {
            dstdata[2] = src1data[0] * src2data[1] - src1data[1] * src2data[0];
            dstdata[0] = src1data[1] * src2data[2] - src1data[2] * src2data[1];
            dstdata[1] = src1data[2] * src2data[0] - src1data[0] * src2data[2];
        }
        else
        {
            int step1 = srcA->step ? srcA->step/sizeof(src1data[0]) : 1;
            int step2 = srcB->step ? srcB->step/sizeof(src1data[0]) : 1;
            int step = dst->step ? dst->step/sizeof(src1data[0]) : 1;

            dstdata[2*step] = src1data[0] * src2data[step2] - src1data[step1] * src2data[0];
            dstdata[0] = src1data[step1] * src2data[step2*2] - src1data[step1*2] * src2data[step2];
            dstdata[step] = src1data[step1*2] * src2data[0] - src1data[0] * src2data[step2*2];
        }
    }
    else if( CV_MAT_DEPTH(type) == CV_64F )
    {
        double* dstdata = (double*)(dst->data.ptr);
        const double* src1data = (double*)(srcA->data.ptr);
        const double* src2data = (double*)(srcB->data.ptr);
        
        if( CV_IS_MAT_CONT(srcA->type & srcB->type & dst->type) )
        {
            dstdata[2] = src1data[0] * src2data[1] - src1data[1] * src2data[0];
            dstdata[0] = src1data[1] * src2data[2] - src1data[2] * src2data[1];
            dstdata[1] = src1data[2] * src2data[0] - src1data[0] * src2data[2];
        }
        else
        {
            int step1 = srcA->step ? srcA->step/sizeof(src1data[0]) : 1;
            int step2 = srcB->step ? srcB->step/sizeof(src1data[0]) : 1;
            int step = dst->step ? dst->step/sizeof(src1data[0]) : 1;

            dstdata[2*step] = src1data[0] * src2data[step2] - src1data[step1] * src2data[0];
            dstdata[0] = src1data[step1] * src2data[step2*2] - src1data[step1*2] * src2data[step2];
            dstdata[step] = src1data[step1*2] * src2data[0] - src1data[0] * src2data[step2*2];
        }
    }
    else
    {
        CV_ERROR( CV_StsUnsupportedFormat, "" );
    }

    __END__;
}


CV_IMPL void
cvCalcPCA( const CvArr* data_arr, CvArr* avg_arr, CvArr* eigenvals, CvArr* eigenvects, int flags )
{
    CvMat* tmp_avg = 0;
    CvMat* tmp_avg_r = 0;
    CvMat* tmp_cov = 0;
    CvMat* tmp_evals = 0;
    CvMat* tmp_evects = 0;
    CvMat* tmp_evects2 = 0;
    CvMat* tmp_data = 0;
    
    CV_FUNCNAME( "cvCalcPCA" );

    __BEGIN__;

    CvMat stub, *data = (CvMat*)data_arr;
    CvMat astub, *avg = (CvMat*)avg_arr;
    CvMat evalstub, *evals = (CvMat*)eigenvals;
    CvMat evectstub, *evects = (CvMat*)eigenvects;
    int covar_flags = CV_COVAR_SCALE;
    int i, len, in_count, count, out_count;

    if( !CV_IS_MAT(data) )
        CV_CALL( data = cvGetMat( data, &stub ));

    if( !CV_IS_MAT(avg) )
        CV_CALL( avg = cvGetMat( avg, &astub ));

    if( !CV_IS_MAT(evals) )
        CV_CALL( evals = cvGetMat( evals, &evalstub ));

    if( !CV_IS_MAT(evects) )
        CV_CALL( evects = cvGetMat( evects, &evectstub ));

    if( CV_MAT_CN(data->type) != 1 || CV_MAT_CN(avg->type) != 1 ||
        CV_MAT_CN(evals->type) != 1 || CV_MAT_CN(evects->type) != 1 )
        CV_ERROR( CV_StsUnsupportedFormat, "All the input and output arrays must be 1-channel" );

    if( CV_MAT_DEPTH(avg->type) < CV_32F || !CV_ARE_DEPTHS_EQ(avg, evals) ||
        !CV_ARE_DEPTHS_EQ(avg, evects) )
        CV_ERROR( CV_StsUnsupportedFormat, "All the output arrays must have the same type, 32fC1 or 64fC1" );

    if( flags & CV_PCA_DATA_AS_COL )
    {
        len = data->rows;
        in_count = data->cols;
        covar_flags |= CV_COVAR_COLS;

        if( avg->cols != 1 || avg->rows != len )
            CV_ERROR( CV_StsBadSize,
            "The mean (average) vector should be data->rows x 1 when CV_PCA_DATA_AS_COL is used" );

        CV_CALL( tmp_avg = cvCreateMat( len, 1, CV_64F ));
    }
    else
    {
        len = data->cols;
        in_count = data->rows;
        covar_flags |= CV_COVAR_ROWS;

        if( avg->rows != 1 || avg->cols != len )
            CV_ERROR( CV_StsBadSize,
            "The mean (average) vector should be 1 x data->cols when CV_PCA_DATA_AS_ROW is used" );

        CV_CALL( tmp_avg = cvCreateMat( 1, len, CV_64F ));
    }

    count = MIN(len, in_count);
    out_count = evals->cols + evals->rows - 1;
    
    if( (evals->cols != 1 && evals->rows != 1) || out_count > count )
        CV_ERROR( CV_StsBadSize,
        "The array of eigenvalues must be 1d vector containing "
        "no more than min(data->rows,data->cols) elements" );

    if( evects->cols != len || evects->rows != out_count )
        CV_ERROR( CV_StsBadSize,
        "The matrix of eigenvalues must have the same number of columns as the input vector length "
        "and the same number of rows as the number of eigenvalues" );

    // "scrambled" way to compute PCA (when cols(A)>rows(A)):
    // B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y
    if( len <= in_count )
        covar_flags |= CV_COVAR_NORMAL;

    if( flags & CV_PCA_USE_AVG ){
        covar_flags |= CV_COVAR_USE_AVG;
		CV_CALL( cvConvert( avg, tmp_avg ) );
	}

    CV_CALL( tmp_cov = cvCreateMat( count, count, CV_64F ));
    CV_CALL( tmp_evals = cvCreateMat( 1, count, CV_64F ));
    CV_CALL( tmp_evects = cvCreateMat( count, count, CV_64F ));

    CV_CALL( cvCalcCovarMatrix( &data_arr, 0, tmp_cov, tmp_avg, covar_flags ));
    CV_CALL( cvSVD( tmp_cov, tmp_evals, tmp_evects, 0, CV_SVD_MODIFY_A + CV_SVD_U_T ));
    tmp_evects->rows = out_count;
    tmp_evals->cols = out_count;
    cvZero( evects );
    cvZero( evals );

    if( covar_flags & CV_COVAR_NORMAL )
    {
        CV_CALL( cvConvert( tmp_evects, evects ));
    }
    else
    {
        // CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A
        // CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A'
        int block_count = 0;

        CV_CALL( tmp_data = cvCreateMat( count, count, CV_64F ));
        CV_CALL( tmp_avg_r = cvCreateMat( count, count, CV_64F ));
        CV_CALL( tmp_evects2 = cvCreateMat( count, count, CV_64F ));

        for( i = 0; i < len; i += block_count )
        {
            CvMat data_part, tdata_part, part, dst_part, avg_part, tmp_avg_part;
            int gemm_flags;

            block_count = MIN( count, len - i );

            if( flags & CV_PCA_DATA_AS_COL )
            {
                cvGetRows( data, &data_part, i, i + block_count );
                cvGetRows( tmp_data, &tdata_part, 0, block_count );
                cvGetRows( tmp_avg, &avg_part, i, i + block_count );
                cvGetRows( tmp_avg_r, &tmp_avg_part, 0, block_count );
                gemm_flags = CV_GEMM_B_T;
            }
            else
            {
                cvGetCols( data, &data_part, i, i + block_count );
                cvGetCols( tmp_data, &tdata_part, 0, block_count );
                cvGetCols( tmp_avg, &avg_part, i, i + block_count );
                cvGetCols( tmp_avg_r, &tmp_avg_part, 0, block_count );
                gemm_flags = 0;
            }

            cvGetCols( tmp_evects2, &part, 0, block_count );
            cvGetCols( evects, &dst_part, i, i + block_count );

            cvConvert( &data_part, &tdata_part );
            cvRepeat( &avg_part, &tmp_avg_part );
            cvSub( &tdata_part, &tmp_avg_part, &tdata_part );
            cvGEMM( tmp_evects, &tdata_part, 1, 0, 0, &part, gemm_flags );
            cvConvert( &part, &dst_part );
        }

        // normalize eigenvectors
        for( i = 0; i < count; i++ )
        {
            CvMat ei;
            cvGetRow( evects, &ei, i );
			cvNormalize( &ei, &ei );
        }
    }

    if( tmp_evals->rows != evals->rows )
        cvReshape( tmp_evals, tmp_evals, 1, evals->rows );
    cvConvert( tmp_evals, evals );
    cvConvert( tmp_avg, avg );

    __END__;

    cvReleaseMat( &tmp_avg );
    cvReleaseMat( &tmp_avg_r );
    cvReleaseMat( &tmp_cov );
    cvReleaseMat( &tmp_evals );
    cvReleaseMat( &tmp_evects );
    cvReleaseMat( &tmp_evects2 );
    cvReleaseMat( &tmp_data );
}


CV_IMPL void
cvProjectPCA( const CvArr* data_arr, const CvArr* avg_arr,
              const CvArr* eigenvects, CvArr* result_arr )
{
    uchar* buffer = 0;
    int local_alloc = 0;
    
    CV_FUNCNAME( "cvProjectPCA" );

    __BEGIN__;

    CvMat stub, *data = (CvMat*)data_arr;
    CvMat astub, *avg = (CvMat*)avg_arr;
    CvMat evectstub, *evects = (CvMat*)eigenvects;
    CvMat rstub, *result = (CvMat*)result_arr;
    CvMat avg_repeated;
    int i, len, in_count;
    int gemm_flags, as_cols, convert_data;
    int block_count0, block_count, buf_size, elem_size;
    uchar* tmp_data_ptr;

    if( !CV_IS_MAT(data) )
        CV_CALL( data = cvGetMat( data, &stub ));

    if( !CV_IS_MAT(avg) )
        CV_CALL( avg = cvGetMat( avg, &astub ));

    if( !CV_IS_MAT(evects) )
        CV_CALL( evects = cvGetMat( evects, &evectstub ));

    if( !CV_IS_MAT(result) )
        CV_CALL( result = cvGetMat( result, &rstub ));

    if( CV_MAT_CN(data->type) != 1 || CV_MAT_CN(avg->type) != 1 )
        CV_ERROR( CV_StsUnsupportedFormat, "All the input and output arrays must be 1-channel" );

    if( CV_MAT_TYPE(avg->type) != CV_32FC1 && CV_MAT_TYPE(avg->type) != CV_64FC1 ||
        !CV_ARE_TYPES_EQ(avg, evects) || !CV_ARE_TYPES_EQ(avg, result) )
        CV_ERROR( CV_StsUnsupportedFormat,
        "All the input and output arrays (except for data) must have the same type, 32fC1 or 64fC1" );

    if( (avg->cols != 1 || avg->rows != data->rows) &&
        (avg->rows != 1 || avg->cols != data->cols) )
        CV_ERROR( CV_StsBadSize,
        "The mean (average) vector should be either 1 x data->cols or data->rows x 1" );

    if( avg->cols == 1 )
    {
        len = data->rows;
        in_count = data->cols;

        gemm_flags = CV_GEMM_A_T + CV_GEMM_B_T;
        as_cols = 1;
    }
    else
    {
        len = data->cols;
        in_count = data->rows;

        gemm_flags = CV_GEMM_B_T;
        as_cols = 0;
    }

    if( evects->cols != len )
        CV_ERROR( CV_StsUnmatchedSizes,
        "Eigenvectors must be stored as rows and be of the same size as input vectors" );

    if( result->cols > evects->rows )
        CV_ERROR( CV_StsOutOfRange,
        "The output matrix of coefficients must have the number of columns "
        "less than or equal to the number of eigenvectors (number of rows in eigenvectors matrix)" );

    evects = cvGetRows( evects, &evectstub, 0, result->cols );

    block_count0 = (1 << 16)/len;
    block_count0 = MAX( block_count0, 4 );
    block_count0 = MIN( block_count0, in_count );
    elem_size = CV_ELEM_SIZE(avg->type);
    convert_data = CV_MAT_DEPTH(data->type) < CV_MAT_DEPTH(avg->type);

    buf_size = block_count0*len*((block_count0 > 1) + 1)*elem_size;

    if( buf_size < CV_MAX_LOCAL_SIZE )
    {
        buffer = (uchar*)cvStackAlloc( buf_size );
        local_alloc = 1;
    }
    else
        CV_CALL( buffer = (uchar*)cvAlloc( buf_size ));

    tmp_data_ptr = buffer;
    if( block_count0 > 1 )
    {
        avg_repeated = cvMat( as_cols ? len : block_count0,
                              as_cols ? block_count0 : len, avg->type, buffer );
        cvRepeat( avg, &avg_repeated );
        tmp_data_ptr += block_count0*len*elem_size;
    }
    else
        avg_repeated = *avg;

    for( i = 0; i < in_count; i += block_count )
    {
        CvMat data_part, norm_data, avg_part, *src = &data_part, out_part;
        
        block_count = MIN( block_count0, in_count - i );
        if( as_