+++ /dev/null
-/*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 */
projects / opencv / blobdiff