1 /*M///////////////////////////////////////////////////////////////////////////////////////
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install,
7 // copy or use the software.
10 // Intel License Agreement
11 // For Open Source Computer Vision Library
13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
14 // Third party copyrights are property of their respective owners.
16 // Redistribution and use in source and binary forms, with or without modification,
17 // are permitted provided that the following conditions are met:
19 // * Redistribution's of source code must retain the above copyright notice,
20 // this list of conditions and the following disclaimer.
22 // * Redistribution's in binary form must reproduce the above copyright notice,
23 // this list of conditions and the following disclaimer in the documentation
24 // and/or other materials provided with the distribution.
26 // * The name of Intel Corporation may not be used to endorse or promote products
27 // derived from this software without specific prior written permission.
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the Intel Corporation or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
44 /****************************************************************************************\
45 * [scaled] Identity matrix initialization *
46 \****************************************************************************************/
49 cvSetIdentity( CvArr* array, CvScalar value )
51 CV_FUNCNAME( "cvSetIdentity" );
55 CvMat stub, *mat = (CvMat*)array;
62 if( !CV_IS_MAT( mat ))
65 CV_CALL( mat = cvGetMat( mat, &stub, &coi ));
67 CV_ERROR( CV_BadCOI, "coi is not supported" );
70 size = cvGetMatSize( mat );
71 len = CV_IMIN( size.width, size.height );
73 type = CV_MAT_TYPE(mat->type);
74 pix_size = CV_ELEM_SIZE(type);
75 size.width *= pix_size;
77 if( CV_IS_MAT_CONT( mat->type ))
79 size.width *= size.height;
87 IPPI_CALL( icvSetZero_8u_C1R( data, step, size ));
90 if( type == CV_32FC1 )
92 float val = (float)value.val[0];
93 float* _data = (float*)data;
94 step /= sizeof(_data[0]);
97 for( i = 0; i < len; i += step )
100 else if( type == CV_64FC1 )
102 double val = value.val[0];
103 double* _data = (double*)data;
104 step /= sizeof(_data[0]);
107 for( i = 0; i < len; i += step )
112 uchar* val_ptr = (uchar*)buf;
113 cvScalarToRawData( &value, buf, type, 0 );
116 for( i = 0; i < len; i += step )
117 for( k = 0; k < pix_size; k++ )
118 data[i+k] = val_ptr[k];
125 /****************************************************************************************\
126 * Trace of the matrix *
127 \****************************************************************************************/
130 cvTrace( const CvArr* array )
132 CvScalar sum = {{0,0,0,0}};
134 CV_FUNCNAME( "cvTrace" );
138 CvMat stub, *mat = 0;
140 if( CV_IS_MAT( array ))
143 int type = CV_MAT_TYPE(mat->type);
144 int size = MIN(mat->rows,mat->cols);
145 uchar* data = mat->data.ptr;
147 if( type == CV_32FC1 )
149 int step = mat->step + sizeof(float);
151 for( ; size--; data += step )
152 sum.val[0] += *(float*)data;
156 if( type == CV_64FC1 )
158 int step = mat->step + sizeof(double);
160 for( ; size--; data += step )
161 sum.val[0] += *(double*)data;
166 CV_CALL( mat = cvGetDiag( array, &stub ));
167 CV_CALL( sum = cvSum( mat ));
175 /****************************************************************************************\
177 \****************************************************************************************/
179 /////////////////// macros for inplace transposition of square matrix ////////////////////
181 #define ICV_DEF_TRANSP_INP_CASE_C1( \
184 arrtype* arr1 = arr; \
185 step /= sizeof(arr[0]); \
189 arr += step, arr1++; \
190 arrtype* arr2 = arr; \
191 arrtype* arr3 = arr1; \
195 arrtype t0 = arr2[0]; \
196 arrtype t1 = arr3[0]; \
203 while( arr2 != arr3 ); \
208 #define ICV_DEF_TRANSP_INP_CASE_C3( \
211 arrtype* arr1 = arr; \
213 step /= sizeof(arr[0]); \
215 for( y = 1; y < len; y++ ) \
217 arr += step, arr1 += 3; \
218 arrtype* arr2 = arr; \
219 arrtype* arr3 = arr1; \
221 for( ; arr2!=arr3; arr2+=3, \
224 arrtype t0 = arr2[0]; \
225 arrtype t1 = arr3[0]; \
241 #define ICV_DEF_TRANSP_INP_CASE_C4( \
244 arrtype* arr1 = arr; \
246 step /= sizeof(arr[0]); \
248 for( y = 1; y < len; y++ ) \
250 arr += step, arr1 += 4; \
251 arrtype* arr2 = arr; \
252 arrtype* arr3 = arr1; \
254 for( ; arr2!=arr3; arr2+=4, \
257 arrtype t0 = arr2[0]; \
258 arrtype t1 = arr3[0]; \
278 //////////////// macros for non-inplace transposition of rectangular matrix //////////////
280 #define ICV_DEF_TRANSP_CASE_C1( arrtype ) \
283 srcstep /= sizeof(src[0]); \
284 dststep /= sizeof(dst[0]); \
286 for( y = 0; y <= size.height - 2; y += 2, \
287 src += 2*srcstep, dst += 2 ) \
289 const arrtype* src1 = src + srcstep; \
290 arrtype* dst1 = dst; \
292 for( x = 0; x <= size.width - 2; \
293 x += 2, dst1 += dststep ) \
295 arrtype t0 = src[x]; \
296 arrtype t1 = src1[x]; \
307 if( x < size.width ) \
309 arrtype t0 = src[x]; \
310 arrtype t1 = src1[x]; \
316 if( y < size.height ) \
318 arrtype* dst1 = dst; \
319 for( x = 0; x <= size.width - 2; \
320 x += 2, dst1 += 2*dststep ) \
322 arrtype t0 = src[x]; \
323 arrtype t1 = src[x + 1]; \
325 dst1[dststep] = t1; \
328 if( x < size.width ) \
330 arrtype t0 = src[x]; \
337 #define ICV_DEF_TRANSP_CASE_C3( arrtype ) \
340 srcstep /= sizeof(src[0]); \
341 dststep /= sizeof(dst[0]); \
343 for( ; size.height--; src+=srcstep, dst+=3 )\
346 arrtype* dst1 = dst; \
348 for( x = 0; x < size.width; x += 3, \
351 arrtype t0 = src[x]; \
352 arrtype t1 = src[x + 1]; \
353 arrtype t2 = src[x + 2]; \
363 #define ICV_DEF_TRANSP_CASE_C4( arrtype ) \
366 srcstep /= sizeof(src[0]); \
367 dststep /= sizeof(dst[0]); \
369 for( ; size.height--; src+=srcstep, dst+=4 )\
372 arrtype* dst1 = dst; \
374 for( x = 0; x < size.width; x += 4, \
377 arrtype t0 = src[x]; \
378 arrtype t1 = src[x + 1]; \
393 #define ICV_DEF_TRANSP_INP_FUNC( flavor, arrtype, cn ) \
394 static CvStatus CV_STDCALL \
395 icvTranspose_##flavor( arrtype* arr, int step, CvSize size )\
397 assert( size.width == size.height ); \
399 ICV_DEF_TRANSP_INP_CASE_C##cn( arrtype, size.width ) \
404 #define ICV_DEF_TRANSP_FUNC( flavor, arrtype, cn ) \
405 static CvStatus CV_STDCALL \
406 icvTranspose_##flavor( const arrtype* src, int srcstep, \
407 arrtype* dst, int dststep, CvSize size )\
409 ICV_DEF_TRANSP_CASE_C##cn( arrtype ) \
414 ICV_DEF_TRANSP_INP_FUNC( 8u_C1IR, uchar, 1 )
415 ICV_DEF_TRANSP_INP_FUNC( 8u_C2IR, ushort, 1 )
416 ICV_DEF_TRANSP_INP_FUNC( 8u_C3IR, uchar, 3 )
417 ICV_DEF_TRANSP_INP_FUNC( 16u_C2IR, int, 1 )
418 ICV_DEF_TRANSP_INP_FUNC( 16u_C3IR, ushort, 3 )
419 ICV_DEF_TRANSP_INP_FUNC( 32s_C2IR, int64, 1 )
420 ICV_DEF_TRANSP_INP_FUNC( 32s_C3IR, int, 3 )
421 ICV_DEF_TRANSP_INP_FUNC( 64s_C2IR, int, 4 )
422 ICV_DEF_TRANSP_INP_FUNC( 64s_C3IR, int64, 3 )
423 ICV_DEF_TRANSP_INP_FUNC( 64s_C4IR, int64, 4 )
426 ICV_DEF_TRANSP_FUNC( 8u_C1R, uchar, 1 )
427 ICV_DEF_TRANSP_FUNC( 8u_C2R, ushort, 1 )
428 ICV_DEF_TRANSP_FUNC( 8u_C3R, uchar, 3 )
429 ICV_DEF_TRANSP_FUNC( 16u_C2R, int, 1 )
430 ICV_DEF_TRANSP_FUNC( 16u_C3R, ushort, 3 )
431 ICV_DEF_TRANSP_FUNC( 32s_C2R, int64, 1 )
432 ICV_DEF_TRANSP_FUNC( 32s_C3R, int, 3 )
433 ICV_DEF_TRANSP_FUNC( 64s_C2R, int, 4 )
434 ICV_DEF_TRANSP_FUNC( 64s_C3R, int64, 3 )
435 ICV_DEF_TRANSP_FUNC( 64s_C4R, int64, 4 )
437 CV_DEF_INIT_PIXSIZE_TAB_2D( Transpose, R )
438 CV_DEF_INIT_PIXSIZE_TAB_2D( Transpose, IR )
441 cvTranspose( const CvArr* srcarr, CvArr* dstarr )
443 static CvBtFuncTable tab, inp_tab;
444 static int inittab = 0;
446 CV_FUNCNAME( "cvTranspose" );
450 CvMat sstub, *src = (CvMat*)srcarr;
451 CvMat dstub, *dst = (CvMat*)dstarr;
457 icvInitTransposeIRTable( &inp_tab );
458 icvInitTransposeRTable( &tab );
462 if( !CV_IS_MAT( src ))
465 CV_CALL( src = cvGetMat( src, &sstub, &coi ));
467 CV_ERROR( CV_BadCOI, "coi is not supported" );
470 type = CV_MAT_TYPE( src->type );
471 pix_size = CV_ELEM_SIZE(type);
472 size = cvGetMatSize( src );
474 if( dstarr == srcarr )
480 if( !CV_IS_MAT( dst ))
483 CV_CALL( dst = cvGetMat( dst, &dstub, &coi ));
486 CV_ERROR( CV_BadCOI, "coi is not supported" );
489 if( !CV_ARE_TYPES_EQ( src, dst ))
490 CV_ERROR( CV_StsUnmatchedFormats, "" );
492 if( size.width != dst->height || size.height != dst->width )
493 CV_ERROR( CV_StsUnmatchedSizes, "" );
496 if( src->data.ptr == dst->data.ptr )
498 if( size.width == size.height )
500 CvFunc2D_1A func = (CvFunc2D_1A)(inp_tab.fn_2d[pix_size]);
503 CV_ERROR( CV_StsUnsupportedFormat, "" );
505 IPPI_CALL( func( src->data.ptr, src->step, size ));
509 if( size.width != 1 && size.height != 1 )
510 CV_ERROR( CV_StsBadSize,
511 "Rectangular matrix can not be transposed inplace" );
513 if( !CV_IS_MAT_CONT( src->type & dst->type ))
514 CV_ERROR( CV_StsBadFlag, "In case of inplace column/row transposition "
515 "both source and destination must be continuous" );
520 CV_SWAP( dst->width, dst->height, t );
521 dst->step = dst->height == 1 ? 0 : pix_size;
527 CvFunc2D_2A func = (CvFunc2D_2A)(tab.fn_2d[pix_size]);
530 CV_ERROR( CV_StsUnsupportedFormat, "" );
532 IPPI_CALL( func( src->data.ptr, src->step,
533 dst->data.ptr, dst->step, size ));
540 /****************************************************************************************\
541 * LU decomposition/back substitution *
542 \****************************************************************************************/
544 #define arrtype float
545 #define temptype double
547 typedef CvStatus (CV_STDCALL * CvLUDecompFunc)( double* A, int stepA, CvSize sizeA,
548 void* B, int stepB, CvSize sizeB,
551 typedef CvStatus (CV_STDCALL * CvLUBackFunc)( double* A, int stepA, CvSize sizeA,
552 void* B, int stepB, CvSize sizeB );
555 #define ICV_DEF_LU_DECOMP_FUNC( flavor, arrtype ) \
556 static CvStatus CV_STDCALL \
557 icvLUDecomp_##flavor( double* A, int stepA, CvSize sizeA, \
558 arrtype* B, int stepB, CvSize sizeB, double* _det ) \
560 int n = sizeA.width; \
564 assert( sizeA.width == sizeA.height ); \
568 assert( sizeA.height == sizeB.height ); \
571 stepA /= sizeof(A[0]); \
572 stepB /= sizeof(B[0]); \
574 for( i = 0; i < n; i++, A += stepA, B += stepB ) \
579 double kval = fabs(A[i]), tval; \
581 /* find the pivot element */ \
582 for( j = i + 1; j < n; j++ ) \
585 tval = fabs(tA[i]); \
603 tA = A + stepA*(k - i); \
606 for( j = i; j < n; j++ ) \
609 CV_SWAP( A[j], tA[j], t ); \
614 tB = B + stepB*(k - i); \
616 for( j = 0; j < m; j++ ) \
619 CV_SWAP( B[j], tB[j], t ); \
628 A[i] = tval; /* to replace division with multiplication in LUBack */ \
630 /* update matrix and the right side of the system */ \
631 for( j = i + 1; j < n; j++ ) \
635 double alpha = -tA[i]*tval; \
637 for( k = i + 1; k < n; k++ ) \
638 tA[k] = tA[k] + alpha*A[k]; \
641 for( k = 0; k < m; k++ ) \
642 tB[k] = (arrtype)(tB[k] + alpha*B[k]); \
653 ICV_DEF_LU_DECOMP_FUNC( 32f, float )
654 ICV_DEF_LU_DECOMP_FUNC( 64f, double )
657 #define ICV_DEF_LU_BACK_FUNC( flavor, arrtype ) \
658 static CvStatus CV_STDCALL \
659 icvLUBack_##flavor( double* A, int stepA, CvSize sizeA, \
660 arrtype* B, int stepB, CvSize sizeB ) \
662 int n = sizeA.width; \
663 int m = sizeB.width, i; \
665 assert( m > 0 && sizeA.width == sizeA.height && \
666 sizeA.height == sizeB.height ); \
667 stepA /= sizeof(A[0]); \
668 stepB /= sizeof(B[0]); \
670 A += stepA*(n - 1); \
671 B += stepB*(n - 1); \
673 for( i = n - 1; i >= 0; i--, A -= stepA ) \
676 for( j = 0; j < m; j++ ) \
678 arrtype* tB = B + j; \
681 for( k = n - 1; k > i; k--, tB -= stepB ) \
684 tB[0] = (arrtype)((tB[0] - x)*A[i]); \
692 ICV_DEF_LU_BACK_FUNC( 32f, float )
693 ICV_DEF_LU_BACK_FUNC( 64f, double )
695 static CvFuncTable lu_decomp_tab, lu_back_tab;
696 static int lu_inittab = 0;
698 static void icvInitLUTable( CvFuncTable* decomp_tab,
699 CvFuncTable* back_tab )
701 decomp_tab->fn_2d[0] = (void*)icvLUDecomp_32f;
702 decomp_tab->fn_2d[1] = (void*)icvLUDecomp_64f;
703 back_tab->fn_2d[0] = (void*)icvLUBack_32f;
704 back_tab->fn_2d[1] = (void*)icvLUBack_64f;
709 /****************************************************************************************\
710 * Determinant of the matrix *
711 \****************************************************************************************/
713 #define det2(m) (m(0,0)*m(1,1) - m(0,1)*m(1,0))
714 #define det3(m) (m(0,0)*(m(1,1)*m(2,2) - m(1,2)*m(2,1)) - \
715 m(0,1)*(m(1,0)*m(2,2) - m(1,2)*m(2,0)) + \
716 m(0,2)*(m(1,0)*m(2,1) - m(1,1)*m(2,0)))
719 cvDet( const CvArr* arr )
725 CV_FUNCNAME( "cvDet" );
729 CvMat stub, *mat = (CvMat*)arr;
732 if( !CV_IS_MAT( mat ))
734 CV_CALL( mat = cvGetMat( mat, &stub ));
737 type = CV_MAT_TYPE( mat->type );
739 if( mat->width != mat->height )
740 CV_ERROR( CV_StsBadSize, "The matrix must be square" );
742 #define Mf( y, x ) ((float*)(m + y*step))[x]
743 #define Md( y, x ) ((double*)(m + y*step))[x]
745 if( mat->width == 2 )
747 uchar* m = mat->data.ptr;
748 int step = mat->step;
750 if( type == CV_32FC1 )
754 else if( type == CV_64FC1 )
760 CV_ERROR( CV_StsUnsupportedFormat, "" );
763 else if( mat->width == 3 )
765 uchar* m = mat->data.ptr;
766 int step = mat->step;
768 if( type == CV_32FC1 )
772 else if( type == CV_64FC1 )
778 CV_ERROR( CV_StsUnsupportedFormat, "" );
781 else if( mat->width == 1 )
783 if( type == CV_32FC1 )
785 result = mat->data.fl[0];
787 else if( type == CV_64FC1 )
789 result = mat->data.db[0];
793 CV_ERROR( CV_StsUnsupportedFormat, "" );
798 CvLUDecompFunc decomp_func;
799 CvSize size = cvGetMatSize( mat );
800 const int worktype = CV_64FC1;
801 int buf_size = size.width*size.height*CV_ELEM_SIZE(worktype);
806 icvInitLUTable( &lu_decomp_tab, &lu_back_tab );
810 if( CV_MAT_CN( type ) != 1 || CV_MAT_DEPTH( type ) < CV_32F )
811 CV_ERROR( CV_StsUnsupportedFormat, "" );
813 if( size.width <= CV_MAX_LOCAL_MAT_SIZE )
815 buffer = (uchar*)cvStackAlloc( buf_size );
820 CV_CALL( buffer = (uchar*)cvAlloc( buf_size ));
823 CV_CALL( cvInitMatHeader( &tmat, size.height, size.width, worktype, buffer ));
824 if( type == worktype )
826 CV_CALL( cvCopy( mat, &tmat ));
829 CV_CALL( cvConvert( mat, &tmat ));
831 decomp_func = (CvLUDecompFunc)(lu_decomp_tab.fn_2d[CV_MAT_DEPTH(worktype)-CV_32F]);
832 assert( decomp_func );
834 IPPI_CALL( decomp_func( tmat.data.db, tmat.step, size, 0, 0, size, &result ));
840 /*icvCheckVector_64f( &result, 1 );*/
844 if( buffer && !local_alloc )
852 /****************************************************************************************\
853 * Inverse (or pseudo-inverse) of the matrix *
854 \****************************************************************************************/
856 #define Sf( y, x ) ((float*)(srcdata + y*srcstep))[x]
857 #define Sd( y, x ) ((double*)(srcdata + y*srcstep))[x]
858 #define Df( y, x ) ((float*)(dstdata + y*dststep))[x]
859 #define Dd( y, x ) ((double*)(dstdata + y*dststep))[x]
862 cvInvert( const CvArr* srcarr, CvArr* dstarr, int method )
872 CV_FUNCNAME( "cvInvert" );
876 CvMat sstub, *src = (CvMat*)srcarr;
877 CvMat dstub, *dst = (CvMat*)dstarr;
880 if( !CV_IS_MAT( src ))
881 CV_CALL( src = cvGetMat( src, &sstub ));
883 if( !CV_IS_MAT( dst ))
884 CV_CALL( dst = cvGetMat( dst, &dstub ));
886 type = CV_MAT_TYPE( src->type );
888 if( method == CV_SVD || method == CV_SVD_SYM )
890 int n = MIN(src->rows,src->cols);
891 if( method == CV_SVD_SYM && src->rows != src->cols )
892 CV_ERROR( CV_StsBadSize, "CV_SVD_SYM method is used for non-square matrix" );
894 CV_CALL( u = cvCreateMat( n, src->rows, src->type ));
895 if( method != CV_SVD_SYM )
896 CV_CALL( v = cvCreateMat( n, src->cols, src->type ));
897 CV_CALL( w = cvCreateMat( n, 1, src->type ));
898 CV_CALL( cvSVD( src, w, u, v, CV_SVD_U_T + CV_SVD_V_T ));
900 if( type == CV_32FC1 )
901 result = w->data.fl[0] >= FLT_EPSILON ?
902 w->data.fl[w->rows-1]/w->data.fl[0] : 0;
904 result = w->data.db[0] >= FLT_EPSILON ?
905 w->data.db[w->rows-1]/w->data.db[0] : 0;
907 CV_CALL( cvSVBkSb( w, u, v ? v : u, 0, dst, CV_SVD_U_T + CV_SVD_V_T ));
910 else if( method != CV_LU )
911 CV_ERROR( CV_StsBadArg, "Unknown inversion method" );
913 if( !CV_ARE_TYPES_EQ( src, dst ))
914 CV_ERROR( CV_StsUnmatchedFormats, "" );
916 if( src->width != src->height )
917 CV_ERROR( CV_StsBadSize, "The matrix must be square" );
919 if( !CV_ARE_SIZES_EQ( src, dst ))
920 CV_ERROR( CV_StsUnmatchedSizes, "" );
922 if( type != CV_32FC1 && type != CV_64FC1 )
923 CV_ERROR( CV_StsUnsupportedFormat, "" );
925 if( src->width <= 3 )
927 uchar* srcdata = src->data.ptr;
928 uchar* dstdata = dst->data.ptr;
929 int srcstep = src->step;
930 int dststep = dst->step;
932 if( src->width == 2 )
934 if( type == CV_32FC1 )
971 else if( src->width == 3 )
973 if( type == CV_32FC1 )
982 t[0] = (float)((Sf(1,1) * Sf(2,2) - Sf(1,2) * Sf(2,1)) * d);
983 t[1] = (float)((Sf(0,2) * Sf(2,1) - Sf(0,1) * Sf(2,2)) * d);
984 t[2] = (float)((Sf(0,1) * Sf(1,2) - Sf(0,2) * Sf(1,1)) * d);
986 t[3] = (float)((Sf(1,2) * Sf(2,0) - Sf(1,0) * Sf(2,2)) * d);
987 t[4] = (float)((Sf(0,0) * Sf(2,2) - Sf(0,2) * Sf(2,0)) * d);
988 t[5] = (float)((Sf(0,2) * Sf(1,0) - Sf(0,0) * Sf(1,2)) * d);
990 t[6] = (float)((Sf(1,0) * Sf(2,1) - Sf(1,1) * Sf(2,0)) * d);
991 t[7] = (float)((Sf(0,1) * Sf(2,0) - Sf(0,0) * Sf(2,1)) * d);
992 t[8] = (float)((Sf(0,0) * Sf(1,1) - Sf(0,1) * Sf(1,0)) * d);
994 Df(0,0) = t[0]; Df(0,1) = t[1]; Df(0,2) = t[2];
995 Df(1,0) = t[3]; Df(1,1) = t[4]; Df(1,2) = t[5];
996 Df(2,0) = t[6]; Df(2,1) = t[7]; Df(2,2) = t[8];
1001 double d = det3(Sd);
1008 t[0] = (Sd(1,1) * Sd(2,2) - Sd(1,2) * Sd(2,1)) * d;
1009 t[1] = (Sd(0,2) * Sd(2,1) - Sd(0,1) * Sd(2,2)) * d;
1010 t[2] = (Sd(0,1) * Sd(1,2) - Sd(0,2) * Sd(1,1)) * d;
1012 t[3] = (Sd(1,2) * Sd(2,0) - Sd(1,0) * Sd(2,2)) * d;
1013 t[4] = (Sd(0,0) * Sd(2,2) - Sd(0,2) * Sd(2,0)) * d;
1014 t[5] = (Sd(0,2) * Sd(1,0) - Sd(0,0) * Sd(1,2)) * d;
1016 t[6] = (Sd(1,0) * Sd(2,1) - Sd(1,1) * Sd(2,0)) * d;
1017 t[7] = (Sd(0,1) * Sd(2,0) - Sd(0,0) * Sd(2,1)) * d;
1018 t[8] = (Sd(0,0) * Sd(1,1) - Sd(0,1) * Sd(1,0)) * d;
1020 Dd(0,0) = t[0]; Dd(0,1) = t[1]; Dd(0,2) = t[2];
1021 Dd(1,0) = t[3]; Dd(1,1) = t[4]; Dd(1,2) = t[5];
1022 Dd(2,0) = t[6]; Dd(2,1) = t[7]; Dd(2,2) = t[8];
1028 assert( src->width == 1 );
1030 if( type == CV_32FC1 )
1036 Df(0,0) = (float)(1./d);
1052 CvLUDecompFunc decomp_func;
1053 CvLUBackFunc back_func;
1054 CvSize size = cvGetMatSize( src );
1055 const int worktype = CV_64FC1;
1056 int buf_size = size.width*size.height*CV_ELEM_SIZE(worktype);
1061 icvInitLUTable( &lu_decomp_tab, &lu_back_tab );
1065 if( size.width <= CV_MAX_LOCAL_MAT_SIZE )
1067 buffer = (uchar*)cvStackAlloc( buf_size );
1072 CV_CALL( buffer = (uchar*)cvAlloc( buf_size ));
1075 CV_CALL( cvInitMatHeader( &tmat, size.height, size.width, worktype, buffer ));
1076 if( type == worktype )
1078 CV_CALL( cvCopy( src, &tmat ));
1081 CV_CALL( cvConvert( src, &tmat ));
1082 CV_CALL( cvSetIdentity( dst ));
1084 decomp_func = (CvLUDecompFunc)(lu_decomp_tab.fn_2d[CV_MAT_DEPTH(type)-CV_32F]);
1085 back_func = (CvLUBackFunc)(lu_back_tab.fn_2d[CV_MAT_DEPTH(type)-CV_32F]);
1086 assert( decomp_func && back_func );
1088 IPPI_CALL( decomp_func( tmat.data.db, tmat.step, size,
1089 dst->data.ptr, dst->step, size, &result ));
1093 IPPI_CALL( back_func( tmat.data.db, tmat.step, size,
1094 dst->data.ptr, dst->step, size ));
1099 CV_CALL( cvSetZero( dst ));
1103 if( buffer && !local_alloc )
1117 /****************************************************************************************\
1118 * Linear system [least-squares] solution *
1119 \****************************************************************************************/
1122 cvSolve( const CvArr* A, const CvArr* b, CvArr* x, int method )
1129 int local_alloc = 0;
1132 CV_FUNCNAME( "cvSolve" );
1136 CvMat sstub, *src = (CvMat*)A;
1137 CvMat dstub, *dst = (CvMat*)x;
1138 CvMat bstub, *src2 = (CvMat*)b;
1141 if( !CV_IS_MAT( src ))
1142 CV_CALL( src = cvGetMat( src, &sstub ));
1144 if( !CV_IS_MAT( src2 ))
1145 CV_CALL( src2 = cvGetMat( src2, &bstub ));
1147 if( !CV_IS_MAT( dst ))
1148 CV_CALL( dst = cvGetMat( dst, &dstub ));
1150 if( method == CV_SVD || method == CV_SVD_SYM )
1152 int n = MIN(src->rows,src->cols);
1154 if( method == CV_SVD_SYM && src->rows != src->cols )
1155 CV_ERROR( CV_StsBadSize, "CV_SVD_SYM method is used for non-square matrix" );
1157 CV_CALL( u = cvCreateMat( n, src->rows, src->type ));
1158 if( method != CV_SVD_SYM )
1159 CV_CALL( v = cvCreateMat( n, src->cols, src->type ));
1160 CV_CALL( w = cvCreateMat( n, 1, src->type ));
1161 CV_CALL( cvSVD( src, w, u, v, CV_SVD_U_T + CV_SVD_V_T ));
1162 CV_CALL( cvSVBkSb( w, u, v ? v : u, src2, dst, CV_SVD_U_T + CV_SVD_V_T ));
1165 else if( method != CV_LU )
1166 CV_ERROR( CV_StsBadArg, "Unknown inversion method" );
1168 type = CV_MAT_TYPE( src->type );
1170 if( !CV_ARE_TYPES_EQ( src, dst ) || !CV_ARE_TYPES_EQ( src, src2 ))
1171 CV_ERROR( CV_StsUnmatchedFormats, "" );
1173 if( src->width != src->height )
1174 CV_ERROR( CV_StsBadSize, "The matrix must be square" );
1176 if( !CV_ARE_SIZES_EQ( src2, dst ) || src->width != src2->height )
1177 CV_ERROR( CV_StsUnmatchedSizes, "" );
1179 if( type != CV_32FC1 && type != CV_64FC1 )
1180 CV_ERROR( CV_StsUnsupportedFormat, "" );
1182 // check case of a single equation and small matrix
1183 if( src->width <= 3 && src2->width == 1 )
1185 #define bf(y) ((float*)(bdata + y*src2step))[0]
1186 #define bd(y) ((double*)(bdata + y*src2step))[0]
1188 uchar* srcdata = src->data.ptr;
1189 uchar* bdata = src2->data.ptr;
1190 uchar* dstdata = dst->data.ptr;
1191 int srcstep = src->step;
1192 int src2step = src2->step;
1193 int dststep = dst->step;
1195 if( src->width == 2 )
1197 if( type == CV_32FC1 )
1199 double d = det2(Sf);
1204 t = (float)((bf(0)*Sf(1,1) - bf(1)*Sf(0,1))*d);
1205 Df(1,0) = (float)((bf(1)*Sf(0,0) - bf(0)*Sf(1,0))*d);
1213 double d = det2(Sd);
1218 t = (bd(0)*Sd(1,1) - bd(1)*Sd(0,1))*d;
1219 Dd(1,0) = (bd(1)*Sd(0,0) - bd(0)*Sd(1,0))*d;
1226 else if( src->width == 3 )
1228 if( type == CV_32FC1 )
1230 double d = det3(Sf);
1237 (bf(0)*(Sf(1,1)*Sf(2,2) - Sf(1,2)*Sf(2,1)) -
1238 Sf(0,1)*(bf(1)*Sf(2,2) - Sf(1,2)*bf(2)) +
1239 Sf(0,2)*(bf(1)*Sf(2,1) - Sf(1,1)*bf(2))));
1242 (Sf(0,0)*(bf(1)*Sf(2,2) - Sf(1,2)*bf(2)) -
1243 bf(0)*(Sf(1,0)*Sf(2,2) - Sf(1,2)*Sf(2,0)) +
1244 Sf(0,2)*(Sf(1,0)*bf(2) - bf(1)*Sf(2,0))));
1247 (Sf(0,0)*(Sf(1,1)*bf(2) - bf(1)*Sf(2,1)) -
1248 Sf(0,1)*(Sf(1,0)*bf(2) - bf(1)*Sf(2,0)) +
1249 bf(0)*(Sf(1,0)*Sf(2,1) - Sf(1,1)*Sf(2,0))));
1260 double d = det3(Sd);
1267 t[0] = ((Sd(1,1) * Sd(2,2) - Sd(1,2) * Sd(2,1))*bd(0) +
1268 (Sd(0,2) * Sd(2,1) - Sd(0,1) * Sd(2,2))*bd(1) +
1269 (Sd(0,1) * Sd(1,2) - Sd(0,2) * Sd(1,1))*bd(2))*d;
1271 t[1] = ((Sd(1,2) * Sd(2,0) - Sd(1,0) * Sd(2,2))*bd(0) +
1272 (Sd(0,0) * Sd(2,2) - Sd(0,2) * Sd(2,0))*bd(1) +
1273 (Sd(0,2) * Sd(1,0) - Sd(0,0) * Sd(1,2))*bd(2))*d;
1275 t[2] = ((Sd(1,0) * Sd(2,1) - Sd(1,1) * Sd(2,0))*bd(0) +
1276 (Sd(0,1) * Sd(2,0) - Sd(0,0) * Sd(2,1))*bd(1) +
1277 (Sd(0,0) * Sd(1,1) - Sd(0,1) * Sd(1,0))*bd(2))*d;
1289 assert( src->width == 1 );
1291 if( type == CV_32FC1 )
1295 Df(0,0) = (float)(bf(0)/d);
1303 Dd(0,0) = (bd(0)/d);
1311 CvLUDecompFunc decomp_func;
1312 CvLUBackFunc back_func;
1313 CvSize size = cvGetMatSize( src );
1314 CvSize dstsize = cvGetMatSize( dst );
1315 int worktype = CV_64FC1;
1316 int buf_size = size.width*size.height*CV_ELEM_SIZE(worktype);
1322 icvInitLUTable( &lu_decomp_tab, &lu_back_tab );
1326 if( size.width <= CV_MAX_LOCAL_MAT_SIZE )
1328 buffer = (uchar*)cvStackAlloc( buf_size );
1333 CV_CALL( buffer = (uchar*)cvAlloc( buf_size ));
1336 CV_CALL( cvInitMatHeader( &tmat, size.height, size.width, worktype, buffer ));
1337 if( type == worktype )
1339 CV_CALL( cvCopy( src, &tmat ));
1342 CV_CALL( cvConvert( src, &tmat ));
1344 if( src2->data.ptr != dst->data.ptr )
1346 CV_CALL( cvCopy( src2, dst ));
1349 decomp_func = (CvLUDecompFunc)(lu_decomp_tab.fn_2d[CV_MAT_DEPTH(type)-CV_32F]);
1350 back_func = (CvLUBackFunc)(lu_back_tab.fn_2d[CV_MAT_DEPTH(type)-CV_32F]);
1351 assert( decomp_func && back_func );
1353 IPPI_CALL( decomp_func( tmat.data.db, tmat.step, size,
1354 dst->data.ptr, dst->step, dstsize, &d ));
1358 IPPI_CALL( back_func( tmat.data.db, tmat.step, size,
1359 dst->data.ptr, dst->step, dstsize ));
1366 CV_CALL( cvSetZero( dst ));
1370 if( buffer && !local_alloc )
1385 /****************************************************************************************\
1386 * 3D vector cross-product *
1387 \****************************************************************************************/
1390 cvCrossProduct( const CvArr* srcAarr, const CvArr* srcBarr, CvArr* dstarr )
1392 CV_FUNCNAME( "cvCrossProduct" );
1396 CvMat stubA, *srcA = (CvMat*)srcAarr;
1397 CvMat stubB, *srcB = (CvMat*)srcBarr;
1398 CvMat dstub, *dst = (CvMat*)dstarr;
1401 if( !CV_IS_MAT(srcA))
1402 CV_CALL( srcA = cvGetMat( srcA, &stubA ));
1404 type = CV_MAT_TYPE( srcA->type );
1406 if( srcA->width*srcA->height*CV_MAT_CN(type) != 3 )
1407 CV_ERROR( CV_StsBadArg, "All the input arrays must be continuous 3-vectors" );
1410 CV_ERROR( CV_StsNullPtr, "" );
1412 if( (srcA->type & ~CV_MAT_CONT_FLAG) == (srcB->type & ~CV_MAT_CONT_FLAG) &&
1413 (srcA->type & ~CV_MAT_CONT_FLAG) == (dst->type & ~CV_MAT_CONT_FLAG) )
1415 if( !srcB->data.ptr || !dst->data.ptr )
1416 CV_ERROR( CV_StsNullPtr, "" );
1420 if( !CV_IS_MAT(srcB))
1421 CV_CALL( srcB = cvGetMat( srcB, &stubB ));
1423 if( !CV_IS_MAT(dst))
1424 CV_CALL( dst = cvGetMat( dst, &dstub ));
1426 if( !CV_ARE_TYPES_EQ( srcA, srcB ) ||
1427 !CV_ARE_TYPES_EQ( srcB, dst ))
1428 CV_ERROR( CV_StsUnmatchedFormats, "" );
1431 if( !CV_ARE_SIZES_EQ( srcA, srcB ) || !CV_ARE_SIZES_EQ( srcB, dst ))
1432 CV_ERROR( CV_StsUnmatchedSizes, "" );
1434 if( CV_MAT_DEPTH(type) == CV_32F )
1436 float* dstdata = (float*)(dst->data.ptr);
1437 const float* src1data = (float*)(srcA->data.ptr);
1438 const float* src2data = (float*)(srcB->data.ptr);
1440 if( CV_IS_MAT_CONT(srcA->type & srcB->type & dst->type) )
1442 dstdata[2] = src1data[0] * src2data[1] - src1data[1] * src2data[0];
1443 dstdata[0] = src1data[1] * src2data[2] - src1data[2] * src2data[1];
1444 dstdata[1] = src1data[2] * src2data[0] - src1data[0] * src2data[2];
1448 int step1 = srcA->step ? srcA->step/sizeof(src1data[0]) : 1;
1449 int step2 = srcB->step ? srcB->step/sizeof(src1data[0]) : 1;
1450 int step = dst->step ? dst->step/sizeof(src1data[0]) : 1;
1452 dstdata[2*step] = src1data[0] * src2data[step2] - src1data[step1] * src2data[0];
1453 dstdata[0] = src1data[step1] * src2data[step2*2] - src1data[step1*2] * src2data[step2];
1454 dstdata[step] = src1data[step1*2] * src2data[0] - src1data[0] * src2data[step2*2];
1457 else if( CV_MAT_DEPTH(type) == CV_64F )
1459 double* dstdata = (double*)(dst->data.ptr);
1460 const double* src1data = (double*)(srcA->data.ptr);
1461 const double* src2data = (double*)(srcB->data.ptr);
1463 if( CV_IS_MAT_CONT(srcA->type & srcB->type & dst->type) )
1465 dstdata[2] = src1data[0] * src2data[1] - src1data[1] * src2data[0];
1466 dstdata[0] = src1data[1] * src2data[2] - src1data[2] * src2data[1];
1467 dstdata[1] = src1data[2] * src2data[0] - src1data[0] * src2data[2];
1471 int step1 = srcA->step ? srcA->step/sizeof(src1data[0]) : 1;
1472 int step2 = srcB->step ? srcB->step/sizeof(src1data[0]) : 1;
1473 int step = dst->step ? dst->step/sizeof(src1data[0]) : 1;
1475 dstdata[2*step] = src1data[0] * src2data[step2] - src1data[step1] * src2data[0];
1476 dstdata[0] = src1data[step1] * src2data[step2*2] - src1data[step1*2] * src2data[step2];
1477 dstdata[step] = src1data[step1*2] * src2data[0] - src1data[0] * src2data[step2*2];
1482 CV_ERROR( CV_StsUnsupportedFormat, "" );
1490 cvCalcPCA( const CvArr* data_arr, CvArr* avg_arr, CvArr* eigenvals, CvArr* eigenvects, int flags )
1493 CvMat* tmp_avg_r = 0;
1495 CvMat* tmp_evals = 0;
1496 CvMat* tmp_evects = 0;
1497 CvMat* tmp_evects2 = 0;
1498 CvMat* tmp_data = 0;
1500 CV_FUNCNAME( "cvCalcPCA" );
1504 CvMat stub, *data = (CvMat*)data_arr;
1505 CvMat astub, *avg = (CvMat*)avg_arr;
1506 CvMat evalstub, *evals = (CvMat*)eigenvals;
1507 CvMat evectstub, *evects = (CvMat*)eigenvects;
1508 int covar_flags = CV_COVAR_SCALE;
1509 int i, len, in_count, count, out_count;
1511 if( !CV_IS_MAT(data) )
1512 CV_CALL( data = cvGetMat( data, &stub ));
1514 if( !CV_IS_MAT(avg) )
1515 CV_CALL( avg = cvGetMat( avg, &astub ));
1517 if( !CV_IS_MAT(evals) )
1518 CV_CALL( evals = cvGetMat( evals, &evalstub ));
1520 if( !CV_IS_MAT(evects) )
1521 CV_CALL( evects = cvGetMat( evects, &evectstub ));
1523 if( CV_MAT_CN(data->type) != 1 || CV_MAT_CN(avg->type) != 1 ||
1524 CV_MAT_CN(evals->type) != 1 || CV_MAT_CN(evects->type) != 1 )
1525 CV_ERROR( CV_StsUnsupportedFormat, "All the input and output arrays must be 1-channel" );
1527 if( CV_MAT_DEPTH(avg->type) < CV_32F || !CV_ARE_DEPTHS_EQ(avg, evals) ||
1528 !CV_ARE_DEPTHS_EQ(avg, evects) )
1529 CV_ERROR( CV_StsUnsupportedFormat, "All the output arrays must have the same type, 32fC1 or 64fC1" );
1531 if( flags & CV_PCA_DATA_AS_COL )
1534 in_count = data->cols;
1535 covar_flags |= CV_COVAR_COLS;
1537 if( avg->cols != 1 || avg->rows != len )
1538 CV_ERROR( CV_StsBadSize,
1539 "The mean (average) vector should be data->rows x 1 when CV_PCA_DATA_AS_COL is used" );
1541 CV_CALL( tmp_avg = cvCreateMat( len, 1, CV_64F ));
1546 in_count = data->rows;
1547 covar_flags |= CV_COVAR_ROWS;
1549 if( avg->rows != 1 || avg->cols != len )
1550 CV_ERROR( CV_StsBadSize,
1551 "The mean (average) vector should be 1 x data->cols when CV_PCA_DATA_AS_ROW is used" );
1553 CV_CALL( tmp_avg = cvCreateMat( 1, len, CV_64F ));
1556 count = MIN(len, in_count);
1557 out_count = evals->cols + evals->rows - 1;
1559 if( (evals->cols != 1 && evals->rows != 1) || out_count > count )
1560 CV_ERROR( CV_StsBadSize,
1561 "The array of eigenvalues must be 1d vector containing "
1562 "no more than min(data->rows,data->cols) elements" );
1564 if( evects->cols != len || evects->rows != out_count )
1565 CV_ERROR( CV_StsBadSize,
1566 "The matrix of eigenvalues must have the same number of columns as the input vector length "
1567 "and the same number of rows as the number of eigenvalues" );
1569 // "scrambled" way to compute PCA (when cols(A)>rows(A)):
1570 // 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
1571 if( len <= in_count )
1572 covar_flags |= CV_COVAR_NORMAL;
1574 if( flags & CV_PCA_USE_AVG ){
1575 covar_flags |= CV_COVAR_USE_AVG;
1576 CV_CALL( cvConvert( avg, tmp_avg ) );
1579 CV_CALL( tmp_cov = cvCreateMat( count, count, CV_64F ));
1580 CV_CALL( tmp_evals = cvCreateMat( 1, count, CV_64F ));
1581 CV_CALL( tmp_evects = cvCreateMat( count, count, CV_64F ));
1583 CV_CALL( cvCalcCovarMatrix( &data_arr, 0, tmp_cov, tmp_avg, covar_flags ));
1584 CV_CALL( cvSVD( tmp_cov, tmp_evals, tmp_evects, 0, CV_SVD_MODIFY_A + CV_SVD_U_T ));
1585 tmp_evects->rows = out_count;
1586 tmp_evals->cols = out_count;
1590 if( covar_flags & CV_COVAR_NORMAL )
1592 CV_CALL( cvConvert( tmp_evects, evects ));
1596 // CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A
1597 // CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A'
1598 int block_count = 0;
1600 CV_CALL( tmp_data = cvCreateMat( count, count, CV_64F ));
1601 CV_CALL( tmp_avg_r = cvCreateMat( count, count, CV_64F ));
1602 CV_CALL( tmp_evects2 = cvCreateMat( out_count, count, CV_64F ));
1604 for( i = 0; i < len; i += block_count )
1606 CvMat data_part, tdata_part, part, dst_part, avg_part, tmp_avg_part;
1609 block_count = MIN( count, len - i );
1611 if( flags & CV_PCA_DATA_AS_COL )
1613 cvGetRows( data, &data_part, i, i + block_count );
1614 cvGetRows( tmp_data, &tdata_part, 0, block_count );
1615 cvGetRows( tmp_avg, &avg_part, i, i + block_count );
1616 cvGetRows( tmp_avg_r, &tmp_avg_part, 0, block_count );
1617 gemm_flags = CV_GEMM_B_T;
1621 cvGetCols( data, &data_part, i, i + block_count );
1622 cvGetCols( tmp_data, &tdata_part, 0, block_count );
1623 cvGetCols( tmp_avg, &avg_part, i, i + block_count );
1624 cvGetCols( tmp_avg_r, &tmp_avg_part, 0, block_count );
1628 cvGetCols( tmp_evects2, &part, 0, block_count );
1629 cvGetCols( evects, &dst_part, i, i + block_count );
1631 cvConvert( &data_part, &tdata_part );
1632 cvRepeat( &avg_part, &tmp_avg_part );
1633 cvSub( &tdata_part, &tmp_avg_part, &tdata_part );
1634 cvGEMM( tmp_evects, &tdata_part, 1, 0, 0, &part, gemm_flags );
1635 cvConvert( &part, &dst_part );
1638 // normalize eigenvectors
1639 for( i = 0; i < out_count; i++ )
1642 cvGetRow( evects, &ei, i );
1643 cvNormalize( &ei, &ei );
1647 if( tmp_evals->rows != evals->rows )
1648 cvReshape( tmp_evals, tmp_evals, 1, evals->rows );
1649 cvConvert( tmp_evals, evals );
1650 cvConvert( tmp_avg, avg );
1654 cvReleaseMat( &tmp_avg );
1655 cvReleaseMat( &tmp_avg_r );
1656 cvReleaseMat( &tmp_cov );
1657 cvReleaseMat( &tmp_evals );
1658 cvReleaseMat( &tmp_evects );
1659 cvReleaseMat( &tmp_evects2 );
1660 cvReleaseMat( &tmp_data );
1665 cvProjectPCA( const CvArr* data_arr, const CvArr* avg_arr,
1666 const CvArr* eigenvects, CvArr* result_arr )
1669 int local_alloc = 0;
1671 CV_FUNCNAME( "cvProjectPCA" );
1675 CvMat stub, *data = (CvMat*)data_arr;
1676 CvMat astub, *avg = (CvMat*)avg_arr;
1677 CvMat evectstub, *evects = (CvMat*)eigenvects;
1678 CvMat rstub, *result = (CvMat*)result_arr;
1680 int i, len, in_count;
1681 int gemm_flags, as_cols, convert_data;
1682 int block_count0, block_count, buf_size, elem_size;
1683 uchar* tmp_data_ptr;
1685 if( !CV_IS_MAT(data) )
1686 CV_CALL( data = cvGetMat( data, &stub ));
1688 if( !CV_IS_MAT(avg) )
1689 CV_CALL( avg = cvGetMat( avg, &astub ));
1691 if( !CV_IS_MAT(evects) )
1692 CV_CALL( evects = cvGetMat( evects, &evectstub ));
1694 if( !CV_IS_MAT(result) )
1695 CV_CALL( result = cvGetMat( result, &rstub ));
1697 if( CV_MAT_CN(data->type) != 1 || CV_MAT_CN(avg->type) != 1 )
1698 CV_ERROR( CV_StsUnsupportedFormat, "All the input and output arrays must be 1-channel" );
1700 if( CV_MAT_TYPE(avg->type) != CV_32FC1 && CV_MAT_TYPE(avg->type) != CV_64FC1 ||
1701 !CV_ARE_TYPES_EQ(avg, evects) || !CV_ARE_TYPES_EQ(avg, result) )
1702 CV_ERROR( CV_StsUnsupportedFormat,
1703 "All the input and output arrays (except for data) must have the same type, 32fC1 or 64fC1" );
1705 if( (avg->cols != 1 || avg->rows != data->rows) &&
1706 (avg->rows != 1 || avg->cols != data->cols) )
1707 CV_ERROR( CV_StsBadSize,
1708 "The mean (average) vector should be either 1 x data->cols or data->rows x 1" );
1710 if( avg->cols == 1 )
1713 in_count = data->cols;
1715 gemm_flags = CV_GEMM_A_T + CV_GEMM_B_T;
1721 in_count = data->rows;
1723 gemm_flags = CV_GEMM_B_T;
1727 if( evects->cols != len )
1728 CV_ERROR( CV_StsUnmatchedSizes,
1729 "Eigenvectors must be stored as rows and be of the same size as input vectors" );
1731 if( result->cols > evects->rows )
1732 CV_ERROR( CV_StsOutOfRange,
1733 "The output matrix of coefficients must have the number of columns "
1734 "less than or equal to the number of eigenvectors (number of rows in eigenvectors matrix)" );
1736 evects = cvGetRows( evects, &evectstub, 0, result->cols );
1738 block_count0 = (1 << 16)/len;
1739 block_count0 = MAX( block_count0, 4 );
1740 block_count0 = MIN( block_count0, in_count );
1741 elem_size = CV_ELEM_SIZE(avg->type);
1742 convert_data = CV_MAT_DEPTH(data->type) < CV_MAT_DEPTH(avg->type);
1744 buf_size = block_count0*len*((block_count0 > 1) + 1)*elem_size;
1746 if( buf_size < CV_MAX_LOCAL_SIZE )
1748 buffer = (uchar*)cvStackAlloc( buf_size );
1752 CV_CALL( buffer = (uchar*)cvAlloc( buf_size ));
1754 tmp_data_ptr = buffer;
1755 if( block_count0 > 1 )
1757 avg_repeated = cvMat( as_cols ? len : block_count0,
1758 as_cols ? block_count0 : len, avg->type, buffer );
1759 cvRepeat( avg, &avg_repeated );
1760 tmp_data_ptr += block_count0*len*elem_size;
1763 avg_repeated = *avg;
1765 for( i = 0; i < in_count; i += block_count )
1767 CvMat data_part, norm_data, avg_part, *src = &data_part, out_part;
1769 block_count = MIN( block_count0, in_count - i );
1772 cvGetCols( data, &data_part, i, i + block_count );
1773 cvGetCols( &avg_repeated, &avg_part, 0, block_count );
1774 norm_data = cvMat( len, block_count, avg->type, tmp_data_ptr );
1778 cvGetRows( data, &data_part, i, i + block_count );
1779 cvGetRows( &avg_repeated, &avg_part, 0, block_count );
1780 norm_data = cvMat( block_count, len, avg->type, tmp_data_ptr );
1785 cvConvert( src, &norm_data );
1789 cvSub( src, &avg_part, &norm_data );
1791 cvGetRows( result, &out_part, i, i + block_count );
1792 cvGEMM( &norm_data, evects, 1, 0, 0, &out_part, gemm_flags );
1803 cvBackProjectPCA( const CvArr* proj_arr, const CvArr* avg_arr,
1804 const CvArr* eigenvects, CvArr* result_arr )
1807 int local_alloc = 0;
1809 CV_FUNCNAME( "cvProjectPCA" );
1813 CvMat pstub, *data = (CvMat*)proj_arr;
1814 CvMat astub, *avg = (CvMat*)avg_arr;
1815 CvMat evectstub, *evects = (CvMat*)eigenvects;
1816 CvMat rstub, *result = (CvMat*)result_arr;
1818 int i, len, in_count, as_cols;
1819 int block_count0, block_count, buf_size, elem_size;
1821 if( !CV_IS_MAT(data) )
1822 CV_CALL( data = cvGetMat( data, &pstub ));
1824 if( !CV_IS_MAT(avg) )
1825 CV_CALL( avg = cvGetMat( avg, &astub ));
1827 if( !CV_IS_MAT(evects) )
1828 CV_CALL( evects = cvGetMat( evects, &evectstub ));
1830 if( !CV_IS_MAT(result) )
1831 CV_CALL( result = cvGetMat( result, &rstub ));
1833 if( CV_MAT_TYPE(avg->type) != CV_32FC1 && CV_MAT_TYPE(avg->type) != CV_64FC1 ||
1834 !CV_ARE_TYPES_EQ(avg, data) || !CV_ARE_TYPES_EQ(avg, evects) || !CV_ARE_TYPES_EQ(avg, result) )
1835 CV_ERROR( CV_StsUnsupportedFormat,
1836 "All the input and output arrays must have the same type, 32fC1 or 64fC1" );
1838 if( (avg->cols != 1 || avg->rows != result->rows) &&
1839 (avg->rows != 1 || avg->cols != result->cols) )
1840 CV_ERROR( CV_StsBadSize,
1841 "The mean (average) vector should be either 1 x result->cols or result->rows x 1" );
1843 if( avg->cols == 1 )
1846 in_count = result->cols;
1852 in_count = result->rows;
1856 if( evects->cols != len )
1857 CV_ERROR( CV_StsUnmatchedSizes,
1858 "Eigenvectors must be stored as rows and be of the same size as the output vectors" );
1860 if( data->cols > evects->rows )
1861 CV_ERROR( CV_StsOutOfRange,
1862 "The input matrix of coefficients must have the number of columns "
1863 "less than or equal to the number of eigenvectors (number of rows in eigenvectors matrix)" );
1865 evects = cvGetRows( evects, &evectstub, 0, data->cols );
1867 block_count0 = (1 << 16)/len;
1868 block_count0 = MAX( block_count0, 4 );
1869 block_count0 = MIN( block_count0, in_count );
1870 elem_size = CV_ELEM_SIZE(avg->type);
1872 buf_size = block_count0*len*(block_count0 > 1)*elem_size;
1874 if( buf_size < CV_MAX_LOCAL_SIZE )
1876 buffer = (uchar*)cvStackAlloc( MAX(buf_size,16) );
1880 CV_CALL( buffer = (uchar*)cvAlloc( buf_size ));
1882 if( block_count0 > 1 )
1884 avg_repeated = cvMat( as_cols ? len : block_count0,
1885 as_cols ? block_count0 : len, avg->type, buffer );
1886 cvRepeat( avg, &avg_repeated );
1889 avg_repeated = *avg;
1891 for( i = 0; i < in_count; i += block_count )
1893 CvMat data_part, avg_part, out_part;
1895 block_count = MIN( block_count0, in_count - i );
1896 cvGetRows( data, &data_part, i, i + block_count );
1900 cvGetCols( result, &out_part, i, i + block_count );
1901 cvGetCols( &avg_repeated, &avg_part, 0, block_count );
1902 cvGEMM( evects, &data_part, 1, &avg_part, 1, &out_part, CV_GEMM_A_T + CV_GEMM_B_T );
1906 cvGetRows( result, &out_part, i, i + block_count );
1907 cvGetRows( &avg_repeated, &avg_part, 0, block_count );
1908 cvGEMM( &data_part, evects, 1, &avg_part, 1, &out_part, 0 );