| /*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*/ |
| |
| #include "_cxcore.h" |
| |
| /*F/////////////////////////////////////////////////////////////////////////////////////// |
| // Names: icvJacobiEigens_32f, icvJacobiEigens_64d |
| // Purpose: Eigenvalues & eigenvectors calculation of a symmetric matrix: |
| // A Vi = Ei Vi |
| // Context: |
| // Parameters: A(n, n) - source symmetric matrix (n - rows & columns number), |
| // V(n, n) - matrix of its eigenvectors |
| // (i-th row is an eigenvector Vi), |
| // E(n) - vector of its eigenvalues |
| // (i-th element is an eigenvalue Ei), |
| // eps - accuracy of diagonalization. |
| // |
| // Returns: |
| // CV_NO_ERROR or error code |
| // Notes: |
| // 1. The functions destroy source matrix A, so if you need it further, you |
| // have to copy it before the processing. |
| // 2. Eigenvalies and eigenvectors are sorted in Ei absolute value descending. |
| // 3. Calculation time depends on eps value. If the time isn't very important, |
| // we recommend to set eps = 0. |
| //F*/ |
| |
| /*=========================== Single precision function ================================*/ |
| |
| static CvStatus CV_STDCALL |
| icvJacobiEigens_32f(float *A, float *V, float *E, int n, float eps) |
| { |
| int i, j, k, ind, iters = 0; |
| float *AA = A, *VV = V; |
| double Amax, anorm = 0, ax; |
| |
| if( A == NULL || V == NULL || E == NULL ) |
| return CV_NULLPTR_ERR; |
| if( n <= 0 ) |
| return CV_BADSIZE_ERR; |
| if( eps < DBL_EPSILON ) |
| eps = DBL_EPSILON; |
| |
| /*-------- Prepare --------*/ |
| for( i = 0; i < n; i++, VV += n, AA += n ) |
| { |
| for( j = 0; j < i; j++ ) |
| { |
| double Am = AA[j]; |
| |
| anorm += Am * Am; |
| } |
| for( j = 0; j < n; j++ ) |
| VV[j] = 0.f; |
| VV[i] = 1.f; |
| } |
| |
| anorm = sqrt( anorm + anorm ); |
| ax = anorm * eps / n; |
| Amax = anorm; |
| |
| while( Amax > ax && iters++ < 100 ) |
| { |
| Amax /= n; |
| do /* while (ind) */ |
| { |
| int p, q; |
| float *V1 = V, *A1 = A; |
| |
| ind = 0; |
| for( p = 0; p < n - 1; p++, A1 += n, V1 += n ) |
| { |
| float *A2 = A + n * (p + 1), *V2 = V + n * (p + 1); |
| |
| for( q = p + 1; q < n; q++, A2 += n, V2 += n ) |
| { |
| double x, y, c, s, c2, s2, a; |
| float *A3, Apq = A1[q], App, Aqq, Aip, Aiq, Vpi, Vqi; |
| |
| if( fabs( Apq ) < Amax ) |
| continue; |
| |
| ind = 1; |
| |
| /*---- Calculation of rotation angle's sine & cosine ----*/ |
| App = A1[p]; |
| Aqq = A2[q]; |
| y = 5.0e-1 * (App - Aqq); |
| x = -Apq / sqrt( (double)Apq * Apq + (double)y * y ); |
| if( y < 0.0 ) |
| x = -x; |
| s = x / sqrt( 2.0 * (1.0 + sqrt( 1.0 - (double)x * x ))); |
| s2 = s * s; |
| c = sqrt( 1.0 - s2 ); |
| c2 = c * c; |
| a = 2.0 * Apq * c * s; |
| |
| /*---- Apq annulation ----*/ |
| A3 = A; |
| for( i = 0; i < p; i++, A3 += n ) |
| { |
| Aip = A3[p]; |
| Aiq = A3[q]; |
| Vpi = V1[i]; |
| Vqi = V2[i]; |
| A3[p] = (float) (Aip * c - Aiq * s); |
| A3[q] = (float) (Aiq * c + Aip * s); |
| V1[i] = (float) (Vpi * c - Vqi * s); |
| V2[i] = (float) (Vqi * c + Vpi * s); |
| } |
| for( ; i < q; i++, A3 += n ) |
| { |
| Aip = A1[i]; |
| Aiq = A3[q]; |
| Vpi = V1[i]; |
| Vqi = V2[i]; |
| A1[i] = (float) (Aip * c - Aiq * s); |
| A3[q] = (float) (Aiq * c + Aip * s); |
| V1[i] = (float) (Vpi * c - Vqi * s); |
| V2[i] = (float) (Vqi * c + Vpi * s); |
| } |
| for( ; i < n; i++ ) |
| { |
| Aip = A1[i]; |
| Aiq = A2[i]; |
| Vpi = V1[i]; |
| Vqi = V2[i]; |
| A1[i] = (float) (Aip * c - Aiq * s); |
| A2[i] = (float) (Aiq * c + Aip * s); |
| V1[i] = (float) (Vpi * c - Vqi * s); |
| V2[i] = (float) (Vqi * c + Vpi * s); |
| } |
| A1[p] = (float) (App * c2 + Aqq * s2 - a); |
| A2[q] = (float) (App * s2 + Aqq * c2 + a); |
| A1[q] = A2[p] = 0.0f; |
| } /*q */ |
| } /*p */ |
| } |
| while( ind ); |
| Amax /= n; |
| } /* while ( Amax > ax ) */ |
| |
| for( i = 0, k = 0; i < n; i++, k += n + 1 ) |
| E[i] = A[k]; |
| /*printf(" M = %d\n", M); */ |
| |
| /* -------- ordering -------- */ |
| for( i = 0; i < n; i++ ) |
| { |
| int m = i; |
| float Em = (float) fabs( E[i] ); |
| |
| for( j = i + 1; j < n; j++ ) |
| { |
| float Ej = (float) fabs( E[j] ); |
| |
| m = (Em < Ej) ? j : m; |
| Em = (Em < Ej) ? Ej : Em; |
| } |
| if( m != i ) |
| { |
| int l; |
| float b = E[i]; |
| |
| E[i] = E[m]; |
| E[m] = b; |
| for( j = 0, k = i * n, l = m * n; j < n; j++, k++, l++ ) |
| { |
| b = V[k]; |
| V[k] = V[l]; |
| V[l] = b; |
| } |
| } |
| } |
| |
| return CV_NO_ERR; |
| } |
| |
| /*=========================== Double precision function ================================*/ |
| |
| static CvStatus CV_STDCALL |
| icvJacobiEigens_64d(double *A, double *V, double *E, int n, double eps) |
| { |
| int i, j, k, p, q, ind, iters = 0; |
| double *A1 = A, *V1 = V, *A2 = A, *V2 = V; |
| double Amax = 0.0, anorm = 0.0, ax; |
| |
| if( A == NULL || V == NULL || E == NULL ) |
| return CV_NULLPTR_ERR; |
| if( n <= 0 ) |
| return CV_BADSIZE_ERR; |
| if( eps < DBL_EPSILON ) |
| eps = DBL_EPSILON; |
| |
| /*-------- Prepare --------*/ |
| for( i = 0; i < n; i++, V1 += n, A1 += n ) |
| { |
| for( j = 0; j < i; j++ ) |
| { |
| double Am = A1[j]; |
| |
| anorm += Am * Am; |
| } |
| for( j = 0; j < n; j++ ) |
| V1[j] = 0.0; |
| V1[i] = 1.0; |
| } |
| |
| anorm = sqrt( anorm + anorm ); |
| ax = anorm * eps / n; |
| Amax = anorm; |
| |
| while( Amax > ax && iters++ < 100 ) |
| { |
| Amax /= n; |
| do /* while (ind) */ |
| { |
| ind = 0; |
| A1 = A; |
| V1 = V; |
| for( p = 0; p < n - 1; p++, A1 += n, V1 += n ) |
| { |
| A2 = A + n * (p + 1); |
| V2 = V + n * (p + 1); |
| for( q = p + 1; q < n; q++, A2 += n, V2 += n ) |
| { |
| double x, y, c, s, c2, s2, a; |
| double *A3, Apq, App, Aqq, App2, Aqq2, Aip, Aiq, Vpi, Vqi; |
| |
| if( fabs( A1[q] ) < Amax ) |
| continue; |
| Apq = A1[q]; |
| |
| ind = 1; |
| |
| /*---- Calculation of rotation angle's sine & cosine ----*/ |
| App = A1[p]; |
| Aqq = A2[q]; |
| y = 5.0e-1 * (App - Aqq); |
| x = -Apq / sqrt( Apq * Apq + (double)y * y ); |
| if( y < 0.0 ) |
| x = -x; |
| s = x / sqrt( 2.0 * (1.0 + sqrt( 1.0 - (double)x * x ))); |
| s2 = s * s; |
| c = sqrt( 1.0 - s2 ); |
| c2 = c * c; |
| a = 2.0 * Apq * c * s; |
| |
| /*---- Apq annulation ----*/ |
| A3 = A; |
| for( i = 0; i < p; i++, A3 += n ) |
| { |
| Aip = A3[p]; |
| Aiq = A3[q]; |
| Vpi = V1[i]; |
| Vqi = V2[i]; |
| A3[p] = Aip * c - Aiq * s; |
| A3[q] = Aiq * c + Aip * s; |
| V1[i] = Vpi * c - Vqi * s; |
| V2[i] = Vqi * c + Vpi * s; |
| } |
| for( ; i < q; i++, A3 += n ) |
| { |
| Aip = A1[i]; |
| Aiq = A3[q]; |
| Vpi = V1[i]; |
| Vqi = V2[i]; |
| A1[i] = Aip * c - Aiq * s; |
| A3[q] = Aiq * c + Aip * s; |
| V1[i] = Vpi * c - Vqi * s; |
| V2[i] = Vqi * c + Vpi * s; |
| } |
| for( ; i < n; i++ ) |
| { |
| Aip = A1[i]; |
| Aiq = A2[i]; |
| Vpi = V1[i]; |
| Vqi = V2[i]; |
| A1[i] = Aip * c - Aiq * s; |
| A2[i] = Aiq * c + Aip * s; |
| V1[i] = Vpi * c - Vqi * s; |
| V2[i] = Vqi * c + Vpi * s; |
| } |
| App2 = App * c2 + Aqq * s2 - a; |
| Aqq2 = App * s2 + Aqq * c2 + a; |
| A1[p] = App2; |
| A2[q] = Aqq2; |
| A1[q] = A2[p] = 0.0; |
| } /*q */ |
| } /*p */ |
| } |
| while( ind ); |
| } /* while ( Amax > ax ) */ |
| |
| for( i = 0, k = 0; i < n; i++, k += n + 1 ) |
| E[i] = A[k]; |
| |
| /* -------- ordering -------- */ |
| for( i = 0; i < n; i++ ) |
| { |
| int m = i; |
| double Em = fabs( E[i] ); |
| |
| for( j = i + 1; j < n; j++ ) |
| { |
| double Ej = fabs( E[j] ); |
| |
| m = (Em < Ej) ? j : m; |
| Em = (Em < Ej) ? Ej : Em; |
| } |
| if( m != i ) |
| { |
| int l; |
| double b = E[i]; |
| |
| E[i] = E[m]; |
| E[m] = b; |
| for( j = 0, k = i * n, l = m * n; j < n; j++, k++, l++ ) |
| { |
| b = V[k]; |
| V[k] = V[l]; |
| V[l] = b; |
| } |
| } |
| } |
| |
| return CV_NO_ERR; |
| } |
| |
| |
| CV_IMPL void |
| cvEigenVV( CvArr* srcarr, CvArr* evectsarr, CvArr* evalsarr, double eps ) |
| { |
| |
| CV_FUNCNAME( "cvEigenVV" ); |
| |
| __BEGIN__; |
| |
| CvMat sstub, *src = (CvMat*)srcarr; |
| CvMat estub1, *evects = (CvMat*)evectsarr; |
| CvMat estub2, *evals = (CvMat*)evalsarr; |
| |
| if( !CV_IS_MAT( src )) |
| CV_CALL( src = cvGetMat( src, &sstub )); |
| |
| if( !CV_IS_MAT( evects )) |
| CV_CALL( evects = cvGetMat( evects, &estub1 )); |
| |
| if( !CV_IS_MAT( evals )) |
| CV_CALL( evals = cvGetMat( evals, &estub2 )); |
| |
| if( src->cols != src->rows ) |
| CV_ERROR( CV_StsUnmatchedSizes, "source is not quadratic matrix" ); |
| |
| if( !CV_ARE_SIZES_EQ( src, evects) ) |
| CV_ERROR( CV_StsUnmatchedSizes, "eigenvectors matrix has inappropriate size" ); |
| |
| if( (evals->rows != src->rows || evals->cols != 1) && |
| (evals->cols != src->rows || evals->rows != 1)) |
| CV_ERROR( CV_StsBadSize, "eigenvalues vector has inappropriate size" ); |
| |
| if( !CV_ARE_TYPES_EQ( src, evects ) || !CV_ARE_TYPES_EQ( src, evals )) |
| CV_ERROR( CV_StsUnmatchedFormats, |
| "input matrix, eigenvalues and eigenvectors must have the same type" ); |
| |
| if( !CV_IS_MAT_CONT( src->type & evals->type & evects->type )) |
| CV_ERROR( CV_BadStep, "all the matrices must be continuous" ); |
| |
| if( CV_MAT_TYPE(src->type) == CV_32FC1 ) |
| { |
| IPPI_CALL( icvJacobiEigens_32f( src->data.fl, |
| evects->data.fl, |
| evals->data.fl, src->cols, (float)eps )); |
| |
| } |
| else if( CV_MAT_TYPE(src->type) == CV_64FC1 ) |
| { |
| IPPI_CALL( icvJacobiEigens_64d( src->data.db, |
| evects->data.db, |
| evals->data.db, src->cols, eps )); |
| } |
| else |
| { |
| CV_ERROR( CV_StsUnsupportedFormat, "Only 32fC1 and 64fC1 types are supported" ); |
| } |
| |
| CV_CHECK_NANS( evects ); |
| CV_CHECK_NANS( evals ); |
| |
| __END__; |
| } |
| |
| /* End of file */ |
