SVD.cpp
Go to the documentation of this file.
1 // ****************************************************************************
2 // This file is part of the Integrating Vision Toolkit (IVT).
3 //
4 // The IVT is maintained by the Karlsruhe Institute of Technology (KIT)
5 // (www.kit.edu) in cooperation with the company Keyetech (www.keyetech.de).
6 //
7 // Copyright (C) 2014 Karlsruhe Institute of Technology (KIT).
8 // All rights reserved.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are met:
12 //
13 // 1. Redistributions of source code must retain the above copyright
14 // notice, this list of conditions and the following disclaimer.
15 //
16 // 2. Redistributions in binary form must reproduce the above copyright
17 // notice, this list of conditions and the following disclaimer in the
18 // documentation and/or other materials provided with the distribution.
19 //
20 // 3. Neither the name of the KIT nor the names of its contributors may be
21 // used to endorse or promote products derived from this software
22 // without specific prior written permission.
23 //
24 // THIS SOFTWARE IS PROVIDED BY THE KIT AND CONTRIBUTORS “AS IS” AND ANY
25 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
26 // WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
27 // DISCLAIMED. IN NO EVENT SHALL THE KIT OR CONTRIBUTORS BE LIABLE FOR ANY
28 // DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
29 // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30 // LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
31 // ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
32 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
33 // THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34 // ****************************************************************************
36 //
37 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
38 //
39 // By downloading, copying, installing or using the software you agree to this license.
40 // If you do not agree to this license, do not download, install,
41 // copy or use the software.
42 //
43 //
44 // Intel License Agreement
45 // For Open Source Computer Vision Library
46 //
47 // Copyright (C) 2000, Intel Corporation, all rights reserved.
48 // Third party copyrights are property of their respective owners.
49 //
50 // Redistribution and use in source and binary forms, with or without modification,
51 // are permitted provided that the following conditions are met:
52 //
53 // * Redistribution's of source code must retain the above copyright notice,
54 // this list of conditions and the following disclaimer.
55 //
56 // * Redistribution's in binary form must reproduce the above copyright notice,
57 // this list of conditions and the following disclaimer in the documentation
58 // and/or other materials provided with the distribution.
59 //
60 // * The name of Intel Corporation may not be used to endorse or promote products
61 // derived from this software without specific prior written permission.
62 //
63 // This software is provided by the copyright holders and contributors "as is" and
64 // any express or implied warranties, including, but not limited to, the implied
65 // warranties of merchantability and fitness for a particular purpose are disclaimed.
66 // In no event shall the Intel Corporation or contributors be liable for any direct,
67 // indirect, incidental, special, exemplary, or consequential damages
68 // (including, but not limited to, procurement of substitute goods or services;
69 // loss of use, data, or profits; or business interruption) however caused
70 // and on any theory of liability, whether in contract, strict liability,
71 // or tort (including negligence or otherwise) arising in any way out of
72 // the use of this software, even if advised of the possibility of such damage.
74 
75 #include <new> // for explicitly using correct new/delete operators on VC DSPs
76 
77 #include "Math/FloatMatrix.h"
78 #include "Math/DoubleMatrix.h"
79 #include "Math/LinearAlgebra.h"
80 #include "Image/ImageProcessor.h"
81 #include "Helpers/helpers.h"
82 
83 #include <math.h>
84 #include <float.h>
85 #include <string.h>
86 #include <stdlib.h>
87 #include <stdio.h>
88 #include <assert.h>
89 #include <algorithm>
90 
91 
93 
94 #define icvGivens_64f( n, x, y, c, s ) \
95 { \
96  int _i; \
97  double* _x = (x); \
98  double* _y = (y); \
99  \
100  for( _i = 0; _i < n; _i++ ) \
101  { \
102  double t0 = _x[_i]; \
103  double t1 = _y[_i]; \
104  _x[_i] = t0*c + t1*s; \
105  _y[_i] = -t0*s + t1*c; \
106  } \
107 }
108 
109 // Taken directly from OpenCV's headers.
110 // CV_CHECK_NANS is just ignored.
111 #define CV_CHECK_NANS( arr )
112 #define CV_SWAP(a,b,t) ((t) = (a), (a) = (b), (b) = (t))
113 
114 // Constants that can be passed as flags to cvSVD.
115 #define CV_SVD_MODIFY_A (1 << 0)
116 #define CV_SVD_U_T (1 << 1)
117 #define CV_SVD_V_T (1 << 2)
118 
119 // This could be #defined to be alloca on UNIX systems.
120 #define cvStackAlloc malloc
121 
122 namespace
123 {
124 
125 /* y[0:m,0:n] += diag(a[0:1,0:m]) * x[0:m,0:n] */
126 void
127 icvMatrAXPY_64f( int m, int n, const double* x, int dx,
128  const double* a, double* y, int dy )
129 {
130  int i, j;
131 
132  for( i = 0; i < m; i++, x += dx, y += dy )
133  {
134  double s = a[i];
135 
136  for( j = 0; j <= n - 4; j += 4 )
137  {
138  double t0 = y[j] + s*x[j];
139  double t1 = y[j+1] + s*x[j+1];
140  y[j] = t0;
141  y[j+1] = t1;
142  t0 = y[j+2] + s*x[j+2];
143  t1 = y[j+3] + s*x[j+3];
144  y[j+2] = t0;
145  y[j+3] = t1;
146  }
147 
148  for( ; j < n; j++ ) y[j] += s*x[j];
149  }
150 }
151 
152 
153 /* y[1:m,-1] = h*y[1:m,0:n]*x[0:1,0:n]'*x[-1] (this is used for U&V reconstruction)
154  y[1:m,0:n] += h*y[1:m,0:n]*x[0:1,0:n]'*x[0:1,0:n] */
155 void
156 icvMatrAXPY3_64f( int m, int n, const double* x, int l, double* y, double h )
157 {
158  int i, j;
159 
160  for( i = 1; i < m; i++ )
161  {
162  double s = 0;
163 
164  y += l;
165 
166  for( j = 0; j <= n - 4; j += 4 )
167  s += x[j]*y[j] + x[j+1]*y[j+1] + x[j+2]*y[j+2] + x[j+3]*y[j+3];
168 
169  for( ; j < n; j++ ) s += x[j]*y[j];
170 
171  s *= h;
172  y[-1] = s*x[-1];
173 
174  for( j = 0; j <= n - 4; j += 4 )
175  {
176  double t0 = y[j] + s*x[j];
177  double t1 = y[j+1] + s*x[j+1];
178  y[j] = t0;
179  y[j+1] = t1;
180  t0 = y[j+2] + s*x[j+2];
181  t1 = y[j+3] + s*x[j+3];
182  y[j+2] = t0;
183  y[j+3] = t1;
184  }
185 
186  for( ; j < n; j++ ) y[j] += s*x[j];
187  }
188 }
189 
190 
191 #define icvGivens_32f( n, x, y, c, s ) \
192 { \
193  int _i; \
194  float* _x = (x); \
195  float* _y = (y); \
196  \
197  for( _i = 0; _i < n; _i++ ) \
198  { \
199  double t0 = _x[_i]; \
200  double t1 = _y[_i]; \
201  _x[_i] = (float)(t0*c + t1*s); \
202  _y[_i] = (float)(-t0*s + t1*c);\
203  } \
204 }
205 
206 void
207 icvMatrAXPY_32f( int m, int n, const float* x, int dx,
208  const float* a, float* y, int dy )
209 {
210  int i, j;
211 
212  for( i = 0; i < m; i++, x += dx, y += dy )
213  {
214  double s = a[i];
215 
216  for( j = 0; j <= n - 4; j += 4 )
217  {
218  double t0 = y[j] + s*x[j];
219  double t1 = y[j+1] + s*x[j+1];
220  y[j] = (float)t0;
221  y[j+1] = (float)t1;
222  t0 = y[j+2] + s*x[j+2];
223  t1 = y[j+3] + s*x[j+3];
224  y[j+2] = (float)t0;
225  y[j+3] = (float)t1;
226  }
227 
228  for( ; j < n; j++ )
229  y[j] = (float)(y[j] + s*x[j]);
230  }
231 }
232 
233 
234 void
235 icvMatrAXPY3_32f( int m, int n, const float* x, int l, float* y, double h )
236 {
237  int i, j;
238 
239  for( i = 1; i < m; i++ )
240  {
241  double s = 0;
242  y += l;
243 
244  for( j = 0; j <= n - 4; j += 4 )
245  s += x[j]*y[j] + x[j+1]*y[j+1] + x[j+2]*y[j+2] + x[j+3]*y[j+3];
246 
247  for( ; j < n; j++ ) s += x[j]*y[j];
248 
249  s *= h;
250  y[-1] = (float)(s*x[-1]);
251 
252  for( j = 0; j <= n - 4; j += 4 )
253  {
254  double t0 = y[j] + s*x[j];
255  double t1 = y[j+1] + s*x[j+1];
256  y[j] = (float)t0;
257  y[j+1] = (float)t1;
258  t0 = y[j+2] + s*x[j+2];
259  t1 = y[j+3] + s*x[j+3];
260  y[j+2] = (float)t0;
261  y[j+3] = (float)t1;
262  }
263 
264  for( ; j < n; j++ ) y[j] = (float)(y[j] + s*x[j]);
265  }
266 }
267 
268 /* accurate hypotenuse calculation */
269 double
270 pythag( double a, double b )
271 {
272  a = fabs( a );
273  b = fabs( b );
274  if( a > b )
275  {
276  b /= a;
277  a *= sqrt( 1. + b * b );
278  }
279  else if( b != 0 )
280  {
281  a /= b;
282  a = b * sqrt( 1. + a * a );
283  }
284 
285  return a;
286 }
287 
288 /****************************************************************************************/
289 /****************************************************************************************/
290 
291 #define MAX_ITERS 30
292 
293 void
294 icvSVD_64f( double* a, int lda, int m, int n,
295  double* w,
296  double* uT, int lduT, int nu,
297  double* vT, int ldvT,
298  double* buffer )
299 {
300  double* e;
301  double* temp;
302  double *w1, *e1;
303  double *hv;
304  double ku0 = 0, kv0 = 0;
305  double anorm = 0;
306  double *a1, *u0 = uT, *v0 = vT;
307  double scale, h;
308  int i, j, k, l;
309  int nm, m1, n1;
310  int nv = n;
311  int iters = 0;
312  double* hv0 = (double*)cvStackAlloc( (m+2)*sizeof(hv0[0])) + 1;
313 
314  e = buffer;
315  w1 = w;
316  e1 = e + 1;
317  nm = n;
318 
319  temp = buffer + nm;
320 
321  memset( w, 0, nm * sizeof( w[0] ));
322  memset( e, 0, nm * sizeof( e[0] ));
323 
324  m1 = m;
325  n1 = n;
326 
327  /* transform a to bi-diagonal form */
328  for( ;; )
329  {
330  int update_u;
331  int update_v;
332 
333  if( m1 == 0 )
334  break;
335 
336  scale = h = 0;
337  update_u = uT && m1 > m - nu;
338  hv = update_u ? uT : hv0;
339 
340  for( j = 0, a1 = a; j < m1; j++, a1 += lda )
341  {
342  double t = a1[0];
343  scale += fabs( hv[j] = t );
344  }
345 
346  if( scale != 0 )
347  {
348  double f = 1./scale, g, s = 0;
349 
350  for( j = 0; j < m1; j++ )
351  {
352  double t = (hv[j] *= f);
353  s += t * t;
354  }
355 
356  g = sqrt( s );
357  f = hv[0];
358  if( f >= 0 )
359  g = -g;
360  hv[0] = f - g;
361  h = 1. / (f * g - s);
362 
363  memset( temp, 0, n1 * sizeof( temp[0] ));
364 
365  /* calc temp[0:n-i] = a[i:m,i:n]'*hv[0:m-i] */
366  icvMatrAXPY_64f( m1, n1 - 1, a + 1, lda, hv, temp + 1, 0 );
367  for( k = 1; k < n1; k++ ) temp[k] *= h;
368 
369  /* modify a: a[i:m,i:n] = a[i:m,i:n] + hv[0:m-i]*temp[0:n-i]' */
370  icvMatrAXPY_64f( m1, n1 - 1, temp + 1, 0, hv, a + 1, lda );
371  *w1 = g*scale;
372  }
373  w1++;
374 
375  /* store -2/(hv'*hv) */
376  if( update_u )
377  {
378  if( m1 == m )
379  ku0 = h;
380  else
381  hv[-1] = h;
382  }
383 
384  a++;
385  n1--;
386  if( vT )
387  vT += ldvT + 1;
388 
389  if( n1 == 0 )
390  break;
391 
392  scale = h = 0;
393  update_v = vT && n1 > n - nv;
394 
395  hv = update_v ? vT : hv0;
396 
397  for( j = 0; j < n1; j++ )
398  {
399  double t = a[j];
400  scale += fabs( hv[j] = t );
401  }
402 
403  if( scale != 0 )
404  {
405  double f = 1./scale, g, s = 0;
406 
407  for( j = 0; j < n1; j++ )
408  {
409  double t = (hv[j] *= f);
410  s += t * t;
411  }
412 
413  g = sqrt( s );
414  f = hv[0];
415  if( f >= 0 )
416  g = -g;
417  hv[0] = f - g;
418  h = 1. / (f * g - s);
419  hv[-1] = 0.;
420 
421  /* update a[i:m:i+1:n] = a[i:m,i+1:n] + (a[i:m,i+1:n]*hv[0:m-i])*... */
422  icvMatrAXPY3_64f( m1, n1, hv, lda, a, h );
423 
424  *e1 = g*scale;
425  }
426  e1++;
427 
428  /* store -2/(hv'*hv) */
429  if( update_v )
430  {
431  if( n1 == n )
432  kv0 = h;
433  else
434  hv[-1] = h;
435  }
436 
437  a += lda;
438  m1--;
439  if( uT )
440  uT += lduT + 1;
441  }
442 
443  m1 -= m1 != 0;
444  n1 -= n1 != 0;
445 
446  /* accumulate left transformations */
447  if( uT )
448  {
449  m1 = m - m1;
450  uT = u0 + m1 * lduT;
451  for( i = m1; i < nu; i++, uT += lduT )
452  {
453  memset( uT + m1, 0, (m - m1) * sizeof( uT[0] ));
454  uT[i] = 1.;
455  }
456 
457  for( i = m1 - 1; i >= 0; i-- )
458  {
459  double s;
460  int lh = nu - i;
461 
462  l = m - i;
463 
464  hv = u0 + (lduT + 1) * i;
465  h = i == 0 ? ku0 : hv[-1];
466 
467  //assert( h <= 0 );
468  if (h > 0) printf("assert: h <= 0 not satisfied, h = %.2f\n", h);
469 
470  if( h != 0 )
471  {
472  uT = hv;
473  icvMatrAXPY3_64f( lh, l-1, hv+1, lduT, uT+1, h );
474 
475  s = hv[0] * h;
476  for( k = 0; k < l; k++ ) hv[k] *= s;
477  hv[0] += 1;
478  }
479  else
480  {
481  for( j = 1; j < l; j++ )
482  hv[j] = 0;
483  for( j = 1; j < lh; j++ )
484  hv[j * lduT] = 0;
485  hv[0] = 1;
486  }
487  }
488  uT = u0;
489  }
490 
491  /* accumulate right transformations */
492  if( vT )
493  {
494  n1 = n - n1;
495  vT = v0 + n1 * ldvT;
496  for( i = n1; i < nv; i++, vT += ldvT )
497  {
498  memset( vT + n1, 0, (n - n1) * sizeof( vT[0] ));
499  vT[i] = 1.;
500  }
501 
502  for( i = n1 - 1; i >= 0; i-- )
503  {
504  double s;
505  int lh = nv - i;
506 
507  l = n - i;
508  hv = v0 + (ldvT + 1) * i;
509  h = i == 0 ? kv0 : hv[-1];
510 
511  //assert( h <= 0 );
512  if (h > 0) printf("assert: h <= 0 not satisfied, h = %.2f (2)\n", h);
513 
514  if( h != 0 )
515  {
516  vT = hv;
517  icvMatrAXPY3_64f( lh, l-1, hv+1, ldvT, vT+1, h );
518 
519  s = hv[0] * h;
520  for( k = 0; k < l; k++ ) hv[k] *= s;
521  hv[0] += 1;
522  }
523  else
524  {
525  for( j = 1; j < l; j++ )
526  hv[j] = 0;
527  for( j = 1; j < lh; j++ )
528  hv[j * ldvT] = 0;
529  hv[0] = 1;
530  }
531  }
532  vT = v0;
533  }
534 
535  for( i = 0; i < nm; i++ )
536  {
537  double tnorm = fabs( w[i] );
538  tnorm += fabs( e[i] );
539 
540  if( anorm < tnorm )
541  anorm = tnorm;
542  }
543 
544  anorm *= DBL_EPSILON;
545 
546  /* diagonalization of the bidiagonal form */
547  for( k = nm - 1; k >= 0; k-- )
548  {
549  double z = 0;
550  iters = 0;
551 
552  for( ;; ) /* do iterations */
553  {
554  double c, s, f, g, x, y;
555  int flag = 0;
556 
557  /* test for splitting */
558  for( l = k; l >= 0; l-- )
559  {
560  if( fabs(e[l]) <= anorm )
561  {
562  flag = 1;
563  break;
564  }
565  //assert( l > 0 );
566  if (l <= 0) printf("assert: l > 0 not satisfied, l = %i\n", l);
567  if( fabs(w[l - 1]) <= anorm )
568  break;
569  }
570 
571  if( !flag )
572  {
573  c = 0;
574  s = 1;
575 
576  for( i = l; i <= k; i++ )
577  {
578  f = s * e[i];
579 
580  e[i] *= c;
581 
582  if( anorm + fabs( f ) == anorm )
583  break;
584 
585  g = w[i];
586  h = pythag( f, g );
587  w[i] = h;
588  c = g / h;
589  s = -f / h;
590 
591  if( uT )
592  icvGivens_64f( m, uT + lduT * (l - 1), uT + lduT * i, c, s );
593  }
594  }
595 
596  z = w[k];
597  if( l == k || iters++ == MAX_ITERS )
598  break;
599 
600  /* shift from bottom 2x2 minor */
601  x = w[l];
602  y = w[k - 1];
603  g = e[k - 1];
604  h = e[k];
605  f = 0.5 * (((g + z) / h) * ((g - z) / y) + y / h - h / y);
606  g = pythag( f, 1 );
607  if( f < 0 )
608  g = -g;
609  f = x - (z / x) * z + (h / x) * (y / (f + g) - h);
610  /* next QR transformation */
611  c = s = 1;
612 
613  for( i = l + 1; i <= k; i++ )
614  {
615  g = e[i];
616  y = w[i];
617  h = s * g;
618  g *= c;
619  z = pythag( f, h );
620  e[i - 1] = z;
621  c = f / z;
622  s = h / z;
623  f = x * c + g * s;
624  g = -x * s + g * c;
625  h = y * s;
626  y *= c;
627 
628  if( vT )
629  icvGivens_64f( n, vT + ldvT * (i - 1), vT + ldvT * i, c, s );
630 
631  z = pythag( f, h );
632  w[i - 1] = z;
633 
634  /* rotation can be arbitrary if z == 0 */
635  if( z != 0 )
636  {
637  c = f / z;
638  s = h / z;
639  }
640  f = c * g + s * y;
641  x = -s * g + c * y;
642 
643  if( uT )
644  icvGivens_64f( m, uT + lduT * (i - 1), uT + lduT * i, c, s );
645  }
646 
647  e[l] = 0;
648  e[k] = f;
649  w[k] = x;
650  } /* end of iteration loop */
651 
652  if( iters > MAX_ITERS )
653  break;
654 
655  if( z < 0 )
656  {
657  w[k] = -z;
658  if( vT )
659  {
660  for( j = 0; j < n; j++ )
661  vT[j + k * ldvT] = -vT[j + k * ldvT];
662  }
663  }
664  } /* end of diagonalization loop */
665 
666  /* sort singular values and corresponding values */
667  for( i = 0; i < nm; i++ )
668  {
669  k = i;
670  for( j = i + 1; j < nm; j++ )
671  if( w[k] < w[j] )
672  k = j;
673 
674  if( k != i )
675  {
676  double t;
677  CV_SWAP( w[i], w[k], t );
678 
679  if( vT )
680  for( j = 0; j < n; j++ )
681  CV_SWAP( vT[j + ldvT*k], vT[j + ldvT*i], t );
682 
683  if( uT )
684  for( j = 0; j < m; j++ )
685  CV_SWAP( uT[j + lduT*k], uT[j + lduT*i], t );
686  }
687  }
688 
689  free(hv0 - 1);
690 }
691 
692 
693 void
694 icvSVD_32f( float* a, int lda, int m, int n,
695  float* w,
696  float* uT, int lduT, int nu,
697  float* vT, int ldvT,
698  float* buffer )
699 {
700  float* e;
701  float* temp;
702  float *w1, *e1;
703  float *hv;
704  double ku0 = 0, kv0 = 0;
705  double anorm = 0;
706  float *a1, *u0 = uT, *v0 = vT;
707  double scale, h;
708  int i, j, k, l;
709  int nm, m1, n1;
710  int nv = n;
711  int iters = 0;
712  float* hv0 = (float*)cvStackAlloc( (m+2)*sizeof(hv0[0])) + 1;
713 
714  e = buffer;
715 
716  w1 = w;
717  e1 = e + 1;
718  nm = n;
719 
720  temp = buffer + nm;
721 
722  memset( w, 0, nm * sizeof( w[0] ));
723  memset( e, 0, nm * sizeof( e[0] ));
724 
725  m1 = m;
726  n1 = n;
727 
728  /* transform a to bi-diagonal form */
729  for( ;; )
730  {
731  int update_u;
732  int update_v;
733 
734  if( m1 == 0 )
735  break;
736 
737  scale = h = 0;
738 
739  update_u = uT && m1 > m - nu;
740  hv = update_u ? uT : hv0;
741 
742  for( j = 0, a1 = a; j < m1; j++, a1 += lda )
743  {
744  double t = a1[0];
745  scale += fabs( hv[j] = (float)t );
746  }
747 
748  if( scale != 0 )
749  {
750  double f = 1./scale, g, s = 0;
751 
752  for( j = 0; j < m1; j++ )
753  {
754  double t = (hv[j] = (float)(hv[j]*f));
755  s += t * t;
756  }
757 
758  g = sqrt( s );
759  f = hv[0];
760  if( f >= 0 )
761  g = -g;
762  hv[0] = (float)(f - g);
763  h = 1. / (f * g - s);
764 
765  memset( temp, 0, n1 * sizeof( temp[0] ));
766 
767  /* calc temp[0:n-i] = a[i:m,i:n]'*hv[0:m-i] */
768  icvMatrAXPY_32f( m1, n1 - 1, a + 1, lda, hv, temp + 1, 0 );
769 
770  for( k = 1; k < n1; k++ ) temp[k] = (float)(temp[k]*h);
771 
772  /* modify a: a[i:m,i:n] = a[i:m,i:n] + hv[0:m-i]*temp[0:n-i]' */
773  icvMatrAXPY_32f( m1, n1 - 1, temp + 1, 0, hv, a + 1, lda );
774  *w1 = (float)(g*scale);
775  }
776  w1++;
777 
778  /* store -2/(hv'*hv) */
779  if( update_u )
780  {
781  if( m1 == m )
782  ku0 = h;
783  else
784  hv[-1] = (float)h;
785  }
786 
787  a++;
788  n1--;
789  if( vT )
790  vT += ldvT + 1;
791 
792  if( n1 == 0 )
793  break;
794 
795  scale = h = 0;
796  update_v = vT && n1 > n - nv;
797  hv = update_v ? vT : hv0;
798 
799  for( j = 0; j < n1; j++ )
800  {
801  double t = a[j];
802  scale += fabs( hv[j] = (float)t );
803  }
804 
805  if( scale != 0 )
806  {
807  double f = 1./scale, g, s = 0;
808 
809  for( j = 0; j < n1; j++ )
810  {
811  double t = (hv[j] = (float)(hv[j]*f));
812  s += t * t;
813  }
814 
815  g = sqrt( s );
816  f = hv[0];
817  if( f >= 0 )
818  g = -g;
819  hv[0] = (float)(f - g);
820  h = 1. / (f * g - s);
821  hv[-1] = 0.f;
822 
823  /* update a[i:m:i+1:n] = a[i:m,i+1:n] + (a[i:m,i+1:n]*hv[0:m-i])*... */
824  icvMatrAXPY3_32f( m1, n1, hv, lda, a, h );
825 
826  *e1 = (float)(g*scale);
827  }
828  e1++;
829 
830  /* store -2/(hv'*hv) */
831  if( update_v )
832  {
833  if( n1 == n )
834  kv0 = h;
835  else
836  hv[-1] = (float)h;
837  }
838 
839  a += lda;
840  m1--;
841  if( uT )
842  uT += lduT + 1;
843  }
844 
845  m1 -= m1 != 0;
846  n1 -= n1 != 0;
847 
848  /* accumulate left transformations */
849  if( uT )
850  {
851  m1 = m - m1;
852  uT = u0 + m1 * lduT;
853  for( i = m1; i < nu; i++, uT += lduT )
854  {
855  memset( uT + m1, 0, (m - m1) * sizeof( uT[0] ));
856  uT[i] = 1.;
857  }
858 
859  for( i = m1 - 1; i >= 0; i-- )
860  {
861  double s;
862  int lh = nu - i;
863 
864  l = m - i;
865 
866  hv = u0 + (lduT + 1) * i;
867  h = i == 0 ? ku0 : hv[-1];
868 
869  //assert( h <= 0 );
870  if (h > 0) printf("assert: h <= 0 not satisfied, h = %.2f (3)\n", h);
871 
872  if( h != 0 )
873  {
874  uT = hv;
875  icvMatrAXPY3_32f( lh, l-1, hv+1, lduT, uT+1, h );
876 
877  s = hv[0] * h;
878  for( k = 0; k < l; k++ ) hv[k] = (float)(hv[k]*s);
879  hv[0] += 1;
880  }
881  else
882  {
883  for( j = 1; j < l; j++ )
884  hv[j] = 0;
885  for( j = 1; j < lh; j++ )
886  hv[j * lduT] = 0;
887  hv[0] = 1;
888  }
889  }
890  uT = u0;
891  }
892 
893  /* accumulate right transformations */
894  if( vT )
895  {
896  n1 = n - n1;
897  vT = v0 + n1 * ldvT;
898  for( i = n1; i < nv; i++, vT += ldvT )
899  {
900  memset( vT + n1, 0, (n - n1) * sizeof( vT[0] ));
901  vT[i] = 1.;
902  }
903 
904  for( i = n1 - 1; i >= 0; i-- )
905  {
906  double s;
907  int lh = nv - i;
908 
909  l = n - i;
910  hv = v0 + (ldvT + 1) * i;
911  h = i == 0 ? kv0 : hv[-1];
912 
913  //assert( h <= 0 );
914  if (h > 0) printf("assert: h <= 0 not satisfied, h = %.2f (4)\n", h);
915 
916  if( h != 0 )
917  {
918  vT = hv;
919  icvMatrAXPY3_32f( lh, l-1, hv+1, ldvT, vT+1, h );
920 
921  s = hv[0] * h;
922  for( k = 0; k < l; k++ ) hv[k] = (float)(hv[k]*s);
923  hv[0] += 1;
924  }
925  else
926  {
927  for( j = 1; j < l; j++ )
928  hv[j] = 0;
929  for( j = 1; j < lh; j++ )
930  hv[j * ldvT] = 0;
931  hv[0] = 1;
932  }
933  }
934  vT = v0;
935  }
936 
937  for( i = 0; i < nm; i++ )
938  {
939  double tnorm = fabs( w[i] );
940  tnorm += fabs( e[i] );
941 
942  if( anorm < tnorm )
943  anorm = tnorm;
944  }
945 
946  anorm *= FLT_EPSILON;
947 
948  /* diagonalization of the bidiagonal form */
949  for( k = nm - 1; k >= 0; k-- )
950  {
951  double z = 0;
952  iters = 0;
953 
954  for( ;; ) /* do iterations */
955  {
956  double c, s, f, g, x, y;
957  int flag = 0;
958 
959  /* test for splitting */
960  for( l = k; l >= 0; l-- )
961  {
962  if( fabs( e[l] ) <= anorm )
963  {
964  flag = 1;
965  break;
966  }
967  //assert( l > 0 );
968  if (l <= 0) printf("assert: l > 0 not satisfied, l = %i\n", l);
969 
970  if( fabs( w[l - 1] ) <= anorm )
971  break;
972  }
973 
974  if( !flag )
975  {
976  c = 0;
977  s = 1;
978 
979  for( i = l; i <= k; i++ )
980  {
981  f = s * e[i];
982  e[i] = (float)(e[i]*c);
983 
984  if( anorm + fabs( f ) == anorm )
985  break;
986 
987  g = w[i];
988  h = pythag( f, g );
989  w[i] = (float)h;
990  c = g / h;
991  s = -f / h;
992 
993  if( uT )
994  icvGivens_32f( m, uT + lduT * (l - 1), uT + lduT * i, c, s );
995  }
996  }
997 
998  z = w[k];
999  if( l == k || iters++ == MAX_ITERS )
1000  break;
1001 
1002  /* shift from bottom 2x2 minor */
1003  x = w[l];
1004  y = w[k - 1];
1005  g = e[k - 1];
1006  h = e[k];
1007  f = 0.5 * (((g + z) / h) * ((g - z) / y) + y / h - h / y);
1008  g = pythag( f, 1 );
1009  if( f < 0 )
1010  g = -g;
1011  f = x - (z / x) * z + (h / x) * (y / (f + g) - h);
1012  /* next QR transformation */
1013  c = s = 1;
1014 
1015  for( i = l + 1; i <= k; i++ )
1016  {
1017  g = e[i];
1018  y = w[i];
1019  h = s * g;
1020  g *= c;
1021  z = pythag( f, h );
1022  e[i - 1] = (float)z;
1023  c = f / z;
1024  s = h / z;
1025  f = x * c + g * s;
1026  g = -x * s + g * c;
1027  h = y * s;
1028  y *= c;
1029 
1030  if( vT )
1031  icvGivens_32f( n, vT + ldvT * (i - 1), vT + ldvT * i, c, s );
1032 
1033  z = pythag( f, h );
1034  w[i - 1] = (float)z;
1035 
1036  /* rotation can be arbitrary if z == 0 */
1037  if( z != 0 )
1038  {
1039  c = f / z;
1040  s = h / z;
1041  }
1042  f = c * g + s * y;
1043  x = -s * g + c * y;
1044 
1045  if( uT )
1046  icvGivens_32f( m, uT + lduT * (i - 1), uT + lduT * i, c, s );
1047  }
1048 
1049  e[l] = 0;
1050  e[k] = (float)f;
1051  w[k] = (float)x;
1052  } /* end of iteration loop */
1053 
1054  if( iters > MAX_ITERS )
1055  break;
1056 
1057  if( z < 0 )
1058  {
1059  w[k] = (float)(-z);
1060  if( vT )
1061  {
1062  for( j = 0; j < n; j++ )
1063  vT[j + k * ldvT] = -vT[j + k * ldvT];
1064  }
1065  }
1066  } /* end of diagonalization loop */
1067 
1068  /* sort singular values and corresponding vectors */
1069  for( i = 0; i < nm; i++ )
1070  {
1071  k = i;
1072  for( j = i + 1; j < nm; j++ )
1073  if( w[k] < w[j] )
1074  k = j;
1075 
1076  if( k != i )
1077  {
1078  float t;
1079  CV_SWAP( w[i], w[k], t );
1080 
1081  if( vT )
1082  for( j = 0; j < n; j++ )
1083  CV_SWAP( vT[j + ldvT*k], vT[j + ldvT*i], t );
1084 
1085  if( uT )
1086  for( j = 0; j < m; j++ )
1087  CV_SWAP( uT[j + lduT*k], uT[j + lduT*i], t );
1088  }
1089  }
1090 
1091  free(hv0 - 1);
1092 }
1093 
1094 
1095 void
1096 icvSVBkSb_64f( int m, int n, const double* w,
1097  const double* uT, int lduT,
1098  const double* vT, int ldvT,
1099  const double* b, int ldb, int nb,
1100  double* x, int ldx, double* buffer )
1101 {
1102  double threshold = 0;
1103  int i, j, nm = MY_MIN( m, n );
1104 
1105  if( !b )
1106  nb = m;
1107 
1108  for( i = 0; i < n; i++ )
1109  memset( x + i*ldx, 0, nb*sizeof(x[0]));
1110 
1111  for( i = 0; i < nm; i++ )
1112  threshold += w[i];
1113  threshold *= 2*DBL_EPSILON;
1114 
1115  /* vT * inv(w) * uT * b */
1116  for( i = 0; i < nm; i++, uT += lduT, vT += ldvT )
1117  {
1118  double wi = w[i];
1119 
1120  if( wi > threshold )
1121  {
1122  wi = 1./wi;
1123 
1124  if( nb == 1 )
1125  {
1126  double s = 0;
1127  if( b )
1128  {
1129  if( ldb == 1 )
1130  {
1131  for( j = 0; j <= m - 4; j += 4 )
1132  s += uT[j]*b[j] + uT[j+1]*b[j+1] + uT[j+2]*b[j+2] + uT[j+3]*b[j+3];
1133  for( ; j < m; j++ )
1134  s += uT[j]*b[j];
1135  }
1136  else
1137  {
1138  for( j = 0; j < m; j++ )
1139  s += uT[j]*b[j*ldb];
1140  }
1141  }
1142  else
1143  s = uT[0];
1144  s *= wi;
1145  if( ldx == 1 )
1146  {
1147  for( j = 0; j <= n - 4; j += 4 )
1148  {
1149  double t0 = x[j] + s*vT[j];
1150  double t1 = x[j+1] + s*vT[j+1];
1151  x[j] = t0;
1152  x[j+1] = t1;
1153  t0 = x[j+2] + s*vT[j+2];
1154  t1 = x[j+3] + s*vT[j+3];
1155  x[j+2] = t0;
1156  x[j+3] = t1;
1157  }
1158 
1159  for( ; j < n; j++ )
1160  x[j] += s*vT[j];
1161  }
1162  else
1163  {
1164  for( j = 0; j < n; j++ )
1165  x[j*ldx] += s*vT[j];
1166  }
1167  }
1168  else
1169  {
1170  if( b )
1171  {
1172  memset( buffer, 0, nb*sizeof(buffer[0]));
1173  icvMatrAXPY_64f( m, nb, b, ldb, uT, buffer, 0 );
1174  for( j = 0; j < nb; j++ )
1175  buffer[j] *= wi;
1176  }
1177  else
1178  {
1179  for( j = 0; j < nb; j++ )
1180  buffer[j] = uT[j]*wi;
1181  }
1182  icvMatrAXPY_64f( n, nb, buffer, 0, vT, x, ldx );
1183  }
1184  }
1185  }
1186 }
1187 
1188 
1189 void
1190 icvSVBkSb_32f( int m, int n, const float* w,
1191  const float* uT, int lduT,
1192  const float* vT, int ldvT,
1193  const float* b, int ldb, int nb,
1194  float* x, int ldx, float* buffer )
1195 {
1196  float threshold = 0.f;
1197  int i, j, nm = MY_MIN( m, n );
1198 
1199  if( !b )
1200  nb = m;
1201 
1202  for( i = 0; i < n; i++ )
1203  memset( x + i*ldx, 0, nb*sizeof(x[0]));
1204 
1205  for( i = 0; i < nm; i++ )
1206  threshold += w[i];
1207  threshold *= 2*FLT_EPSILON;
1208 
1209  /* vT * inv(w) * uT * b */
1210  for( i = 0; i < nm; i++, uT += lduT, vT += ldvT )
1211  {
1212  double wi = w[i];
1213 
1214  if( wi > threshold )
1215  {
1216  wi = 1./wi;
1217 
1218  if( nb == 1 )
1219  {
1220  double s = 0;
1221  if( b )
1222  {
1223  if( ldb == 1 )
1224  {
1225  for( j = 0; j <= m - 4; j += 4 )
1226  s += uT[j]*b[j] + uT[j+1]*b[j+1] + uT[j+2]*b[j+2] + uT[j+3]*b[j+3];
1227  for( ; j < m; j++ )
1228  s += uT[j]*b[j];
1229  }
1230  else
1231  {
1232  for( j = 0; j < m; j++ )
1233  s += uT[j]*b[j*ldb];
1234  }
1235  }
1236  else
1237  s = uT[0];
1238  s *= wi;
1239 
1240  if( ldx == 1 )
1241  {
1242  for( j = 0; j <= n - 4; j += 4 )
1243  {
1244  double t0 = x[j] + s*vT[j];
1245  double t1 = x[j+1] + s*vT[j+1];
1246  x[j] = (float)t0;
1247  x[j+1] = (float)t1;
1248  t0 = x[j+2] + s*vT[j+2];
1249  t1 = x[j+3] + s*vT[j+3];
1250  x[j+2] = (float)t0;
1251  x[j+3] = (float)t1;
1252  }
1253 
1254  for( ; j < n; j++ )
1255  x[j] = (float)(x[j] + s*vT[j]);
1256  }
1257  else
1258  {
1259  for( j = 0; j < n; j++ )
1260  x[j*ldx] = (float)(x[j*ldx] + s*vT[j]);
1261  }
1262  }
1263  else
1264  {
1265  if( b )
1266  {
1267  memset( buffer, 0, nb*sizeof(buffer[0]));
1268  icvMatrAXPY_32f( m, nb, b, ldb, uT, buffer, 0 );
1269  for( j = 0; j < nb; j++ )
1270  buffer[j] = (float)(buffer[j]*wi);
1271  }
1272  else
1273  {
1274  for( j = 0; j < nb; j++ )
1275  buffer[j] = (float)(uT[j]*wi);
1276  }
1277  icvMatrAXPY_32f( n, nb, buffer, 0, vT, x, ldx );
1278  }
1279  }
1280  }
1281 }
1282 
1283 
1284 void
1285 cvSVD( const CFloatMatrix* aarr, CFloatMatrix* warr, CFloatMatrix* uarr, CFloatMatrix* varr, int flags )
1286 {
1287  typedef unsigned char uchar;
1288 
1289  uchar* buffer = 0;
1290  //int local_alloc = 0;
1291 
1292  CFloatMatrix tmat;
1293  CFloatMatrix ustub, vstub;
1294 
1295  const CFloatMatrix *a = aarr;
1296  CFloatMatrix *w = warr;
1297 
1298  uchar* tw = 0;
1299  int a_buf_offset = 0, u_buf_offset = 0, buf_size, pix_size;
1300  int temp_u = 0, /* temporary storage for U is needed */
1301  t_svd; /* special case: a->rows < a->columns */
1302  int m, n;
1303  int w_rows, w_cols;
1304  int u_rows = 0, u_cols = 0;
1305  int w_is_mat = 0;
1306 
1307  if( a->rows >= a->columns )
1308  {
1309  m = a->rows;
1310  n = a->columns;
1311  w_rows = w->rows;
1312  w_cols = w->columns;
1313  t_svd = 0;
1314  }
1315  else
1316  {
1317  CFloatMatrix* t;
1318  CV_SWAP( uarr, varr, t );
1319 
1320  flags = (flags & CV_SVD_U_T ? CV_SVD_V_T : 0)|
1321  (flags & CV_SVD_V_T ? CV_SVD_U_T : 0);
1322  m = a->columns;
1323  n = a->rows;
1324  w_rows = w->columns;
1325  w_cols = w->rows;
1326  t_svd = 1;
1327  }
1328 
1329  CFloatMatrix *u = uarr;
1330  CFloatMatrix *v = varr;
1331 
1332  w_is_mat = w_cols > 1 && w_rows > 1;
1333  if (!w_is_mat && w_cols + w_rows - 1 == n)
1334  tw = (uchar*) w->data;
1335 
1336  if( u )
1337  {
1338  if( !(flags & CV_SVD_U_T) )
1339  {
1340  u_rows = u->rows;
1341  u_cols = u->columns;
1342  }
1343  else
1344  {
1345  u_rows = u->columns;
1346  u_cols = u->rows;
1347  }
1348 
1349  if( u_rows != m || (u_cols != m && u_cols != n))
1350  !t_svd ? printf( "U matrix has unappropriate size" ) : printf( "V matrix has unappropriate size" );
1351 
1352  temp_u = (u_rows != u_cols && !(flags & CV_SVD_U_T)) || u->data == a->data;
1353 
1354  if( w_is_mat && u_cols != w_rows )
1355  !t_svd ? printf( "U and W have incompatible sizes" ) : printf( "V and W have incompatible sizes" );
1356  }
1357  else
1358  {
1359  u = &ustub;
1360  u->rows = 0;
1361  u->columns = 0;
1362  u->data = 0;
1363  }
1364 
1365  if( v )
1366  {
1367  int v_rows, v_cols;
1368 
1369  if( !(flags & CV_SVD_V_T) )
1370  {
1371  v_rows = v->rows;
1372  v_cols = v->columns;
1373  }
1374  else
1375  {
1376  v_rows = v->columns;
1377  v_cols = v->rows;
1378  }
1379 
1380  if( v_rows != n || v_cols != n )
1381  t_svd ? printf( "U matrix has unappropriate size") : printf("V matrix has unappropriate size" );
1382 
1383  if( w_is_mat && w_cols != v_cols )
1384  t_svd ? printf( "U and W have incompatible sizes") : printf("V and W have incompatible sizes" );
1385  }
1386  else
1387  {
1388  v = &vstub;
1389  v->rows = 0;
1390  v->columns = 0;
1391  v->data = 0;
1392  }
1393 
1394  pix_size = sizeof(float);
1395  buf_size = n*2 + m;
1396 
1397  if( !(flags & CV_SVD_MODIFY_A) )
1398  {
1399  a_buf_offset = buf_size;
1400  buf_size += a->rows*a->columns;
1401  }
1402 
1403  if( temp_u )
1404  {
1405  u_buf_offset = buf_size;
1406  buf_size += u->rows*u->columns;
1407  }
1408 
1409  buf_size *= pix_size;
1410 
1411  buffer = (uchar*)malloc( buf_size );
1412 
1413  if( !(flags & CV_SVD_MODIFY_A) )
1414  {
1415  tmat.rows = m;
1416  tmat.columns = n;
1417  tmat.data = (float *)(buffer + a_buf_offset * pix_size);
1418 
1419  if( !t_svd )
1420  ImageProcessor::CopyMatrix(a, &tmat); // cvCopy( a, &tmat );
1421  else
1422  LinearAlgebra::Transpose(a, &tmat); // cvT( a, &tmat );
1423  a = &tmat;
1424  }
1425 
1426  if( temp_u )
1427  {
1428  ustub.rows = u_cols;
1429  ustub.columns = u_rows;
1430  ustub.data = (float *)(buffer + u_buf_offset * pix_size);
1431  u = &ustub;
1432  }
1433 
1434  if( !tw )
1435  tw = buffer + (n + m)*pix_size;
1436 
1437 
1438 
1439  icvSVD_32f( a->data, a->columns, a->rows, a->columns,
1440  (float*)tw, u->data, u->columns, u_cols,
1441  v->data, v->columns, (float*)buffer );
1442 
1443 
1444  if( (void *) tw != (void *) w->data )
1445  {
1446  int shift = w->columns != 1;
1447  ImageProcessor::Zero(w); // cvSetZero( w );
1448 
1449  for( int i = 0; i < n; i++ )
1450  ((float*)(((unsigned char *) w->data) + i*(w->columns*sizeof(float))))[i*shift] = ((float*)tw)[i];
1451  }
1452 
1453  if( uarr )
1454  {
1455  if( !(flags & CV_SVD_U_T))
1456  LinearAlgebra::Transpose(u, uarr); // cvT( u, uarr );
1457  else if( temp_u )
1458  ImageProcessor::CopyMatrix(u, uarr); // cvCopy( u, uarr );
1459  }
1460 
1461  if( varr )
1462  {
1463  if( !(flags & CV_SVD_V_T))
1464  LinearAlgebra::Transpose(v, varr); // cvT( v, varr );
1465  }
1466 
1467  if( buffer )
1468  free( buffer );
1469 }
1470 
1471 void
1472 cvSVD( const CDoubleMatrix* aarr, CDoubleMatrix* warr, CDoubleMatrix* uarr, CDoubleMatrix* varr, int flags )
1473 {
1474  typedef unsigned char uchar;
1475 
1476  uchar* buffer = 0;
1477  //int local_alloc = 0;
1478 
1479  CDoubleMatrix tmat;
1480  CDoubleMatrix ustub, vstub;
1481 
1482  const CDoubleMatrix *a = aarr;
1483  CDoubleMatrix *w = warr;
1484 
1485  uchar* tw = 0;
1486  int a_buf_offset = 0, u_buf_offset = 0, buf_size, pix_size;
1487  int temp_u = 0, /* temporary storage for U is needed */
1488  t_svd; /* special case: a->rows < a->columns */
1489  int m, n;
1490  int w_rows, w_cols;
1491  int u_rows = 0, u_cols = 0;
1492  int w_is_mat = 0;
1493 
1494  if( a->rows >= a->columns )
1495  {
1496  m = a->rows;
1497  n = a->columns;
1498  w_rows = w->rows;
1499  w_cols = w->columns;
1500  t_svd = 0;
1501  }
1502  else
1503  {
1504  CDoubleMatrix* t;
1505  CV_SWAP( uarr, varr, t );
1506 
1507  flags = (flags & CV_SVD_U_T ? CV_SVD_V_T : 0)|
1508  (flags & CV_SVD_V_T ? CV_SVD_U_T : 0);
1509  m = a->columns;
1510  n = a->rows;
1511  w_rows = w->columns;
1512  w_cols = w->rows;
1513  t_svd = 1;
1514  }
1515 
1516  CDoubleMatrix *u = uarr;
1517  CDoubleMatrix *v = varr;
1518 
1519  w_is_mat = w_cols > 1 && w_rows > 1;
1520  if (!w_is_mat && w_cols + w_rows - 1 == n)
1521  tw = (uchar*) w->data;
1522 
1523  if( u )
1524  {
1525  if( !(flags & CV_SVD_U_T) )
1526  {
1527  u_rows = u->rows;
1528  u_cols = u->columns;
1529  }
1530  else
1531  {
1532  u_rows = u->columns;
1533  u_cols = u->rows;
1534  }
1535 
1536  if( u_rows != m || (u_cols != m && u_cols != n))
1537  !t_svd ? printf( "U matrix has unappropriate size" ) : printf( "V matrix has unappropriate size" );
1538 
1539  temp_u = (u_rows != u_cols && !(flags & CV_SVD_U_T)) || u->data == a->data;
1540 
1541  if( w_is_mat && u_cols != w_rows )
1542  !t_svd ? printf( "U and W have incompatible sizes" ) : printf( "V and W have incompatible sizes" );
1543  }
1544  else
1545  {
1546  u = &ustub;
1547  u->rows = 0;
1548  u->columns = 0;
1549  u->data = 0;
1550  }
1551 
1552  if( v )
1553  {
1554  int v_rows, v_cols;
1555 
1556  if( !(flags & CV_SVD_V_T) )
1557  {
1558  v_rows = v->rows;
1559  v_cols = v->columns;
1560  }
1561  else
1562  {
1563  v_rows = v->columns;
1564  v_cols = v->rows;
1565  }
1566 
1567  if( v_rows != n || v_cols != n )
1568  t_svd ? printf( "U matrix has unappropriate size" ) : printf( "V matrix has unappropriate size" );
1569 
1570  if( w_is_mat && w_cols != v_cols )
1571  t_svd ? printf( "U and W have incompatible sizes" ) : printf( "V and W have incompatible sizes" );
1572  }
1573  else
1574  {
1575  v = &vstub;
1576  v->rows = 0;
1577  v->columns = 0;
1578  v->data = 0;
1579  }
1580 
1581  pix_size = sizeof(double);
1582  buf_size = n*2 + m;
1583 
1584  if( !(flags & CV_SVD_MODIFY_A) )
1585  {
1586  a_buf_offset = buf_size;
1587  buf_size += a->rows*a->columns;
1588  }
1589 
1590  if( temp_u )
1591  {
1592  u_buf_offset = buf_size;
1593  buf_size += u->rows*u->columns;
1594  }
1595 
1596  buf_size *= pix_size;
1597 
1598  buffer = (uchar*)malloc( buf_size );
1599 
1600  if( !(flags & CV_SVD_MODIFY_A) )
1601  {
1602  tmat.rows = m;
1603  tmat.columns = n;
1604  tmat.data = (double *)(buffer + a_buf_offset * pix_size);
1605 
1606  if( !t_svd )
1607  ImageProcessor::CopyMatrix(a, &tmat); // cvCopy( a, &tmat );
1608  else
1609  LinearAlgebra::Transpose(a, &tmat); // cvT( a, &tmat );
1610  a = &tmat;
1611  }
1612 
1613  if( temp_u )
1614  {
1615  ustub.rows = u_cols;
1616  ustub.columns = u_rows;
1617  ustub.data = (double *)(buffer + u_buf_offset * pix_size);
1618  u = &ustub;
1619  }
1620 
1621  if( !tw )
1622  tw = buffer + (n + m)*pix_size;
1623 
1624 
1625 
1626  icvSVD_64f( a->data, a->columns, a->rows, a->columns,
1627  (double*)tw, u->data, u->columns, u_cols,
1628  v->data, v->columns, (double*)buffer );
1629 
1630 
1631  if( (void *) tw != (void *) w->data )
1632  {
1633  int shift = w->columns != 1;
1634  ImageProcessor::Zero(w); // cvSetZero( w );
1635 
1636  for( int i = 0; i < n; i++ )
1637  ((double*)(((unsigned char *) w->data) + i*(w->columns*sizeof(double))))[i*shift] = ((double*)tw)[i];
1638  }
1639 
1640  if( uarr )
1641  {
1642  if( !(flags & CV_SVD_U_T))
1643  LinearAlgebra::Transpose(u, uarr); // cvT( u, uarr );
1644  else if( temp_u )
1645  ImageProcessor::CopyMatrix(u, uarr); // cvCopy( u, uarr );
1646  }
1647 
1648  if( varr )
1649  {
1650  if( !(flags & CV_SVD_V_T))
1651  LinearAlgebra::Transpose(v, varr); // cvT( v, varr );
1652  }
1653 
1654  if( buffer )
1655  free( buffer );
1656 }
1657 } // End of unnamed namespace
1658 
1659 
1660 void LinearAlgebra::SVD(const CFloatMatrix *A, CFloatMatrix *W, CFloatMatrix *U, CFloatMatrix *V, bool bAllowModifyA, bool bReturnUTransposed, bool bReturnVTransposed)
1661 {
1662  const int columns = A->columns;
1663  const int rows = A->rows;
1664 
1665  if (W->columns != columns || W->rows != rows)
1666  {
1667  printf("error: W should have %i columns and %i rows for LinearAlgebra::SVD\n", columns, rows);
1668  return;
1669  }
1670 
1671  if (U && (U->columns != rows || U->rows != rows))
1672  {
1673  printf("error: U should have %i columns and %i rows for LinearAlgebra::SVD\n", rows, rows);
1674  return;
1675  }
1676 
1677  if (V && (V->columns != columns || V->rows != columns))
1678  {
1679  printf("error: V should have %i columns and %i rows for LinearAlgebra::SVD\n", columns, columns);
1680  return;
1681  }
1682 
1683  int flags = 0;
1684 
1685  if (bAllowModifyA)
1686  flags |= CV_SVD_MODIFY_A;
1687 
1688  if (bReturnUTransposed)
1689  flags |= CV_SVD_U_T;
1690 
1691  if (bReturnVTransposed)
1692  flags |= CV_SVD_V_T;
1693 
1694  cvSVD(A, W, U, V, flags);
1695 }
1696 
1697 void LinearAlgebra::SVD(const CDoubleMatrix *A, CDoubleMatrix *W, CDoubleMatrix *U, CDoubleMatrix *V, bool bAllowModifyA, bool bReturnUTransposed, bool bReturnVTransposed)
1698 {
1699  const int columns = A->columns;
1700  const int rows = A->rows;
1701 
1702  if (W->columns != columns || W->rows != rows)
1703  {
1704  printf("error: W should have %i columns and %i rows for LinearAlgebra::SVD\n", columns, rows);
1705  return;
1706  }
1707 
1708  if (U && (U->columns != rows || U->rows != rows))
1709  {
1710  printf("error: U should have %i columns and %i rows for LinearAlgebra::SVD\n", rows, rows);
1711  return;
1712  }
1713 
1714  if (V && (V->columns != columns || V->rows != columns))
1715  {
1716  printf("error: V should have %i columns and %i rows for LinearAlgebra::SVD\n", columns, columns);
1717  return;
1718  }
1719 
1720  int flags = 0;
1721 
1722  if (bAllowModifyA)
1723  flags |= CV_SVD_MODIFY_A;
1724 
1725  if (bReturnUTransposed)
1726  flags |= CV_SVD_U_T;
1727 
1728  if (bReturnVTransposed)
1729  flags |= CV_SVD_V_T;
1730 
1731  cvSVD(A, W, U, V, flags);
1732 }
g
GLubyte g
Definition: glext.h:5166
t
GLdouble GLdouble t
Definition: glext.h:3219
CDoubleMatrix::data
double * data
Definition: DoubleMatrix.h:85
z
GLdouble GLdouble z
Definition: glext.h:3343
scale
GLenum GLenum GLenum GLenum GLenum scale
Definition: glext.h:5308
n
GLenum GLsizei n
Definition: glext.h:4209
s
GLdouble s
Definition: glext.h:3211
v0
GLfloat v0
Definition: glext.h:3530
cvStackAlloc
#define cvStackAlloc
Definition: SVD.cpp:120
CV_SVD_V_T
#define CV_SVD_V_T
Definition: SVD.cpp:117
MAX_ITERS
#define MAX_ITERS
Definition: SVD.cpp:291
LinearAlgebra::SVD
void SVD(const CFloatMatrix *A, CFloatMatrix *W, CFloatMatrix *U=0, CFloatMatrix *V=0, bool bAllowModifyA=false, bool bReturnUTransposed=false, bool bReturnVTransposed=false)
Definition: SVD.cpp:1660
x
GLenum GLint x
Definition: glext.h:3125
icvGivens_32f
#define icvGivens_32f(n, x, y, c, s)
Definition: SVD.cpp:191
c
const GLubyte * c
Definition: glext.h:5181
CV_SVD_MODIFY_A
#define CV_SVD_MODIFY_A
Definition: SVD.cpp:115
ImageProcessor::CopyMatrix
bool CopyMatrix(const CFloatMatrix *pInputImage, CFloatMatrix *pOutputImage)
Copies one CFloatMatrix to another.
Definition: ImageProcessor.cpp:2912
CFloatMatrix
Data structure for the representation of a matrix of values of the data type float.
Definition: FloatMatrix.h:56
CV_SVD_U_T
#define CV_SVD_U_T
Definition: SVD.cpp:116
CFloatMatrix::data
float * data
Definition: FloatMatrix.h:91
a
GLubyte GLubyte GLubyte a
Definition: glext.h:5166
CV_SWAP
#define CV_SWAP(a, b, t)
Definition: SVD.cpp:112
b
GLubyte GLubyte b
Definition: glext.h:5166
MY_MIN
#define MY_MIN(a, b)
Definition: helpers.h:64
y
GLenum GLint GLint y
Definition: glext.h:3125
v
const GLdouble * v
Definition: glext.h:3212
CDoubleMatrix
Data structure for the representation of a matrix of values of the data type double.
Definition: DoubleMatrix.h:54
LinearAlgebra::Transpose
void Transpose(const CFloatMatrix *pMatrix, CFloatMatrix *pResultMatrix)
Definition: LinearAlgebra.cpp:254
icvGivens_64f
#define icvGivens_64f(n, x, y, c, s)
Definition: SVD.cpp:94
buffer
GLuint buffer
Definition: glext.h:3384
ImageProcessor::Zero
void Zero(CByteImage *pImage, const MyRegion *pROI=0)
Sets all values in a CByteImage to zero.
Definition: ImageProcessor.cpp:2939
w
GLubyte GLubyte GLubyte GLubyte w
Definition: glext.h:3571
w1
GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint GLdouble w1
Definition: glext.h:4584

asr_ivt
Author(s): Allgeyer Tobias, Hutmacher Robin, Kleinert Daniel, Meißner Pascal, Scholz Jonas, Stöckle Patrick
autogenerated on Mon Dec 2 2019 03:47:28