1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
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.
8 //
9 //
10 // Intel License Agreement
11 // For Open Source Computer Vision Library
12 //
13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
14 // Third party copyrights are property of their respective owners.
15 //
16 // Redistribution and use in source and binary forms, with or without modification,
17 // are permitted provided that the following conditions are met:
18 //
19 // * Redistribution's of source code must retain the above copyright notice,
20 // this list of conditions and the following disclaimer.
21 //
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.
25 //
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.
28 //
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.
39 //
40 //M*/
41
42 #include "_cxcore.h"
43
44 /*F///////////////////////////////////////////////////////////////////////////////////////
45 // Names: icvJacobiEigens_32f, icvJacobiEigens_64d
46 // Purpose: Eigenvalues & eigenvectors calculation of a symmetric matrix:
47 // A Vi = Ei Vi
48 // Context:
49 // Parameters: A(n, n) - source symmetric matrix (n - rows & columns number),
50 // V(n, n) - matrix of its eigenvectors
51 // (i-th row is an eigenvector Vi),
52 // E(n) - vector of its eigenvalues
53 // (i-th element is an eigenvalue Ei),
54 // eps - accuracy of diagonalization.
55 //
56 // Returns:
57 // CV_NO_ERROR or error code
58 // Notes:
59 // 1. The functions destroy source matrix A, so if you need it further, you
60 // have to copy it before the processing.
61 // 2. Eigenvalies and eigenvectors are sorted in Ei absolute value descending.
62 // 3. Calculation time depends on eps value. If the time isn't very important,
63 // we recommend to set eps = 0.
64 //F*/
65
66 /*=========================== Single precision function ================================*/
67
68 static CvStatus CV_STDCALL
icvJacobiEigens_32f(float * A,float * V,float * E,int n,float eps)69 icvJacobiEigens_32f(float *A, float *V, float *E, int n, float eps)
70 {
71 int i, j, k, ind, iters = 0;
72 float *AA = A, *VV = V;
73 double Amax, anorm = 0, ax;
74
75 if( A == NULL || V == NULL || E == NULL )
76 return CV_NULLPTR_ERR;
77 if( n <= 0 )
78 return CV_BADSIZE_ERR;
79 if( eps < DBL_EPSILON )
80 eps = DBL_EPSILON;
81
82 /*-------- Prepare --------*/
83 for( i = 0; i < n; i++, VV += n, AA += n )
84 {
85 for( j = 0; j < i; j++ )
86 {
87 double Am = AA[j];
88
89 anorm += Am * Am;
90 }
91 for( j = 0; j < n; j++ )
92 VV[j] = 0.f;
93 VV[i] = 1.f;
94 }
95
96 anorm = sqrt( anorm + anorm );
97 ax = anorm * eps / n;
98 Amax = anorm;
99
100 while( Amax > ax && iters++ < 100 )
101 {
102 Amax /= n;
103 do /* while (ind) */
104 {
105 int p, q;
106 float *V1 = V, *A1 = A;
107
108 ind = 0;
109 for( p = 0; p < n - 1; p++, A1 += n, V1 += n )
110 {
111 float *A2 = A + n * (p + 1), *V2 = V + n * (p + 1);
112
113 for( q = p + 1; q < n; q++, A2 += n, V2 += n )
114 {
115 double x, y, c, s, c2, s2, a;
116 float *A3, Apq = A1[q], App, Aqq, Aip, Aiq, Vpi, Vqi;
117
118 if( fabs( Apq ) < Amax )
119 continue;
120
121 ind = 1;
122
123 /*---- Calculation of rotation angle's sine & cosine ----*/
124 App = A1[p];
125 Aqq = A2[q];
126 y = 5.0e-1 * (App - Aqq);
127 x = -Apq / sqrt( (double)Apq * Apq + (double)y * y );
128 if( y < 0.0 )
129 x = -x;
130 s = x / sqrt( 2.0 * (1.0 + sqrt( 1.0 - (double)x * x )));
131 s2 = s * s;
132 c = sqrt( 1.0 - s2 );
133 c2 = c * c;
134 a = 2.0 * Apq * c * s;
135
136 /*---- Apq annulation ----*/
137 A3 = A;
138 for( i = 0; i < p; i++, A3 += n )
139 {
140 Aip = A3[p];
141 Aiq = A3[q];
142 Vpi = V1[i];
143 Vqi = V2[i];
144 A3[p] = (float) (Aip * c - Aiq * s);
145 A3[q] = (float) (Aiq * c + Aip * s);
146 V1[i] = (float) (Vpi * c - Vqi * s);
147 V2[i] = (float) (Vqi * c + Vpi * s);
148 }
149 for( ; i < q; i++, A3 += n )
150 {
151 Aip = A1[i];
152 Aiq = A3[q];
153 Vpi = V1[i];
154 Vqi = V2[i];
155 A1[i] = (float) (Aip * c - Aiq * s);
156 A3[q] = (float) (Aiq * c + Aip * s);
157 V1[i] = (float) (Vpi * c - Vqi * s);
158 V2[i] = (float) (Vqi * c + Vpi * s);
159 }
160 for( ; i < n; i++ )
161 {
162 Aip = A1[i];
163 Aiq = A2[i];
164 Vpi = V1[i];
165 Vqi = V2[i];
166 A1[i] = (float) (Aip * c - Aiq * s);
167 A2[i] = (float) (Aiq * c + Aip * s);
168 V1[i] = (float) (Vpi * c - Vqi * s);
169 V2[i] = (float) (Vqi * c + Vpi * s);
170 }
171 A1[p] = (float) (App * c2 + Aqq * s2 - a);
172 A2[q] = (float) (App * s2 + Aqq * c2 + a);
173 A1[q] = A2[p] = 0.0f;
174 } /*q */
175 } /*p */
176 }
177 while( ind );
178 Amax /= n;
179 } /* while ( Amax > ax ) */
180
181 for( i = 0, k = 0; i < n; i++, k += n + 1 )
182 E[i] = A[k];
183 /*printf(" M = %d\n", M); */
184
185 /* -------- ordering -------- */
186 for( i = 0; i < n; i++ )
187 {
188 int m = i;
189 float Em = (float) fabs( E[i] );
190
191 for( j = i + 1; j < n; j++ )
192 {
193 float Ej = (float) fabs( E[j] );
194
195 m = (Em < Ej) ? j : m;
196 Em = (Em < Ej) ? Ej : Em;
197 }
198 if( m != i )
199 {
200 int l;
201 float b = E[i];
202
203 E[i] = E[m];
204 E[m] = b;
205 for( j = 0, k = i * n, l = m * n; j < n; j++, k++, l++ )
206 {
207 b = V[k];
208 V[k] = V[l];
209 V[l] = b;
210 }
211 }
212 }
213
214 return CV_NO_ERR;
215 }
216
217 /*=========================== Double precision function ================================*/
218
219 static CvStatus CV_STDCALL
icvJacobiEigens_64d(double * A,double * V,double * E,int n,double eps)220 icvJacobiEigens_64d(double *A, double *V, double *E, int n, double eps)
221 {
222 int i, j, k, p, q, ind, iters = 0;
223 double *A1 = A, *V1 = V, *A2 = A, *V2 = V;
224 double Amax = 0.0, anorm = 0.0, ax;
225
226 if( A == NULL || V == NULL || E == NULL )
227 return CV_NULLPTR_ERR;
228 if( n <= 0 )
229 return CV_BADSIZE_ERR;
230 if( eps < DBL_EPSILON )
231 eps = DBL_EPSILON;
232
233 /*-------- Prepare --------*/
234 for( i = 0; i < n; i++, V1 += n, A1 += n )
235 {
236 for( j = 0; j < i; j++ )
237 {
238 double Am = A1[j];
239
240 anorm += Am * Am;
241 }
242 for( j = 0; j < n; j++ )
243 V1[j] = 0.0;
244 V1[i] = 1.0;
245 }
246
247 anorm = sqrt( anorm + anorm );
248 ax = anorm * eps / n;
249 Amax = anorm;
250
251 while( Amax > ax && iters++ < 100 )
252 {
253 Amax /= n;
254 do /* while (ind) */
255 {
256 ind = 0;
257 A1 = A;
258 V1 = V;
259 for( p = 0; p < n - 1; p++, A1 += n, V1 += n )
260 {
261 A2 = A + n * (p + 1);
262 V2 = V + n * (p + 1);
263 for( q = p + 1; q < n; q++, A2 += n, V2 += n )
264 {
265 double x, y, c, s, c2, s2, a;
266 double *A3, Apq, App, Aqq, App2, Aqq2, Aip, Aiq, Vpi, Vqi;
267
268 if( fabs( A1[q] ) < Amax )
269 continue;
270 Apq = A1[q];
271
272 ind = 1;
273
274 /*---- Calculation of rotation angle's sine & cosine ----*/
275 App = A1[p];
276 Aqq = A2[q];
277 y = 5.0e-1 * (App - Aqq);
278 x = -Apq / sqrt( Apq * Apq + (double)y * y );
279 if( y < 0.0 )
280 x = -x;
281 s = x / sqrt( 2.0 * (1.0 + sqrt( 1.0 - (double)x * x )));
282 s2 = s * s;
283 c = sqrt( 1.0 - s2 );
284 c2 = c * c;
285 a = 2.0 * Apq * c * s;
286
287 /*---- Apq annulation ----*/
288 A3 = A;
289 for( i = 0; i < p; i++, A3 += n )
290 {
291 Aip = A3[p];
292 Aiq = A3[q];
293 Vpi = V1[i];
294 Vqi = V2[i];
295 A3[p] = Aip * c - Aiq * s;
296 A3[q] = Aiq * c + Aip * s;
297 V1[i] = Vpi * c - Vqi * s;
298 V2[i] = Vqi * c + Vpi * s;
299 }
300 for( ; i < q; i++, A3 += n )
301 {
302 Aip = A1[i];
303 Aiq = A3[q];
304 Vpi = V1[i];
305 Vqi = V2[i];
306 A1[i] = Aip * c - Aiq * s;
307 A3[q] = Aiq * c + Aip * s;
308 V1[i] = Vpi * c - Vqi * s;
309 V2[i] = Vqi * c + Vpi * s;
310 }
311 for( ; i < n; i++ )
312 {
313 Aip = A1[i];
314 Aiq = A2[i];
315 Vpi = V1[i];
316 Vqi = V2[i];
317 A1[i] = Aip * c - Aiq * s;
318 A2[i] = Aiq * c + Aip * s;
319 V1[i] = Vpi * c - Vqi * s;
320 V2[i] = Vqi * c + Vpi * s;
321 }
322 App2 = App * c2 + Aqq * s2 - a;
323 Aqq2 = App * s2 + Aqq * c2 + a;
324 A1[p] = App2;
325 A2[q] = Aqq2;
326 A1[q] = A2[p] = 0.0;
327 } /*q */
328 } /*p */
329 }
330 while( ind );
331 } /* while ( Amax > ax ) */
332
333 for( i = 0, k = 0; i < n; i++, k += n + 1 )
334 E[i] = A[k];
335
336 /* -------- ordering -------- */
337 for( i = 0; i < n; i++ )
338 {
339 int m = i;
340 double Em = fabs( E[i] );
341
342 for( j = i + 1; j < n; j++ )
343 {
344 double Ej = fabs( E[j] );
345
346 m = (Em < Ej) ? j : m;
347 Em = (Em < Ej) ? Ej : Em;
348 }
349 if( m != i )
350 {
351 int l;
352 double b = E[i];
353
354 E[i] = E[m];
355 E[m] = b;
356 for( j = 0, k = i * n, l = m * n; j < n; j++, k++, l++ )
357 {
358 b = V[k];
359 V[k] = V[l];
360 V[l] = b;
361 }
362 }
363 }
364
365 return CV_NO_ERR;
366 }
367
368
369 CV_IMPL void
cvEigenVV(CvArr * srcarr,CvArr * evectsarr,CvArr * evalsarr,double eps)370 cvEigenVV( CvArr* srcarr, CvArr* evectsarr, CvArr* evalsarr, double eps )
371 {
372
373 CV_FUNCNAME( "cvEigenVV" );
374
375 __BEGIN__;
376
377 CvMat sstub, *src = (CvMat*)srcarr;
378 CvMat estub1, *evects = (CvMat*)evectsarr;
379 CvMat estub2, *evals = (CvMat*)evalsarr;
380
381 if( !CV_IS_MAT( src ))
382 CV_CALL( src = cvGetMat( src, &sstub ));
383
384 if( !CV_IS_MAT( evects ))
385 CV_CALL( evects = cvGetMat( evects, &estub1 ));
386
387 if( !CV_IS_MAT( evals ))
388 CV_CALL( evals = cvGetMat( evals, &estub2 ));
389
390 if( src->cols != src->rows )
391 CV_ERROR( CV_StsUnmatchedSizes, "source is not quadratic matrix" );
392
393 if( !CV_ARE_SIZES_EQ( src, evects) )
394 CV_ERROR( CV_StsUnmatchedSizes, "eigenvectors matrix has inappropriate size" );
395
396 if( (evals->rows != src->rows || evals->cols != 1) &&
397 (evals->cols != src->rows || evals->rows != 1))
398 CV_ERROR( CV_StsBadSize, "eigenvalues vector has inappropriate size" );
399
400 if( !CV_ARE_TYPES_EQ( src, evects ) || !CV_ARE_TYPES_EQ( src, evals ))
401 CV_ERROR( CV_StsUnmatchedFormats,
402 "input matrix, eigenvalues and eigenvectors must have the same type" );
403
404 if( !CV_IS_MAT_CONT( src->type & evals->type & evects->type ))
405 CV_ERROR( CV_BadStep, "all the matrices must be continuous" );
406
407 if( CV_MAT_TYPE(src->type) == CV_32FC1 )
408 {
409 IPPI_CALL( icvJacobiEigens_32f( src->data.fl,
410 evects->data.fl,
411 evals->data.fl, src->cols, (float)eps ));
412
413 }
414 else if( CV_MAT_TYPE(src->type) == CV_64FC1 )
415 {
416 IPPI_CALL( icvJacobiEigens_64d( src->data.db,
417 evects->data.db,
418 evals->data.db, src->cols, eps ));
419 }
420 else
421 {
422 CV_ERROR( CV_StsUnsupportedFormat, "Only 32fC1 and 64fC1 types are supported" );
423 }
424
425 CV_CHECK_NANS( evects );
426 CV_CHECK_NANS( evals );
427
428 __END__;
429 }
430
431 /* End of file */
432