fixed bug in SVD, ticket #2027; fixed building highgui with ffmpeg support on MacOSX
authorVadim Pisarevsky <no@email>
Thu, 28 Jun 2012 19:45:13 +0000 (19:45 +0000)
committerVadim Pisarevsky <no@email>
Thu, 28 Jun 2012 19:45:13 +0000 (19:45 +0000)
modules/core/src/lapack.cpp
modules/core/test/test_operations.cpp
modules/highgui/CMakeLists.txt

index 12ad278..c66923d 100644 (file)
@@ -527,7 +527,7 @@ template<> inline int VBLAS<double>::givensx(double* a, double* b, int n, double
 #endif
 
 template<typename _Tp> void
-JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, int m, int n, int n1)
+JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, int m, int n, int n1, double minval)
 {
     VBLAS<_Tp> vblas;
     AutoBuffer<double> Wbuf(n);
@@ -587,11 +587,11 @@ JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, int m, int
                     delta = p*p*0.5/(gamma + beta);
                 }
                 
-                if( iter % 2 )
+                W[i] += delta;
+                W[j] -= delta;
+                
+                if( iter % 2 != 0 && W[i] > 0 && W[j] > 0 )
                 {
-                    W[i] += delta;
-                    W[j] -= delta;
-                    
                     k = vblas.givens(Ai, Aj, m, c, s);
                     
                     for( ; k < m; k++ )
@@ -671,12 +671,13 @@ JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, int m, int
     
     if( !Vt )
         return;
+    
     RNG rng(0x12345678);
     for( i = 0; i < n1; i++ )
     {
         sd = i < n ? W[i] : 0;
         
-        while( sd == 0 )
+        while( sd <= minval )
         {
             // if we got a zero singular value, then in order to get the corresponding left singular vector
             // we generate a random vector, project it to the previously computed left singular vectors,
@@ -724,12 +725,12 @@ JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, int m, int
     
 static void JacobiSVD(float* At, size_t astep, float* W, float* Vt, size_t vstep, int m, int n, int n1=-1)
 {
-    JacobiSVDImpl_(At, astep, W, Vt, vstep, m, n, !Vt ? 0 : n1 < 0 ? n : n1);
+    JacobiSVDImpl_(At, astep, W, Vt, vstep, m, n, !Vt ? 0 : n1 < 0 ? n : n1, FLT_MIN);
 }
 
 static void JacobiSVD(double* At, size_t astep, double* W, double* Vt, size_t vstep, int m, int n, int n1=-1)
 {
-    JacobiSVDImpl_(At, astep, W, Vt, vstep, m, n, !Vt ? 0 : n1 < 0 ? n : n1);
+    JacobiSVDImpl_(At, astep, W, Vt, vstep, m, n, !Vt ? 0 : n1 < 0 ? n : n1, DBL_MIN);
 }
     
 /* y[0:m,0:n] += diag(a[0:1,0:m]) * x[0:m,0:n] */
index fd9c42b..04c6098 100644 (file)
@@ -949,6 +949,34 @@ bool CV_OperationsTest::TestSVD()
         SVD svd(A, SVD::FULL_UV); 
         if( norm(A*svd.vt.row(3).t(), CV_C) > FLT_EPSILON )
             throw test_excep();
+        
+        Mat Dp(3,3,CV_32FC1);
+        Mat Dc(3,3,CV_32FC1);
+        Mat Q(3,3,CV_32FC1);
+        Mat U,Vt,R,T,W;
+        
+        Dp.at<float>(0,0)=0.86483884; Dp.at<float>(0,1)= -0.3077251; Dp.at<float>(0,2)=-0.55711365;
+        Dp.at<float>(1,0)=0.49294353; Dp.at<float>(1,1)=-0.24209651; Dp.at<float>(1,2)=-0.25084701;
+        Dp.at<float>(2,0)=0;          Dp.at<float>(2,1)=0;           Dp.at<float>(2,2)=0;
+        
+        Dc.at<float>(0,0)=0.75632739; Dc.at<float>(0,1)= -0.38859656; Dc.at<float>(0,2)=-0.36773083;
+        Dc.at<float>(1,0)=0.9699229; Dc.at<float>(1,1)=-0.49858192; Dc.at<float>(1,2)=-0.47134098;
+        Dc.at<float>(2,0)=0.10566688; Dc.at<float>(2,1)=-0.060333252; Dc.at<float>(2,2)=-0.045333147;
+        
+        Q=Dp*Dc.t();
+        SVD decomp;
+        decomp=SVD(Q);
+        U=decomp.u;
+        Vt=decomp.vt;
+        W=decomp.w;
+        Mat I = Mat::eye(3, 3, CV_32F);
+        
+        if( norm(U*U.t(), I, CV_C) > FLT_EPSILON ||
+            norm(Vt*Vt.t(), I, CV_C) > FLT_EPSILON ||
+            W.at<float>(2) < 0 || W.at<float>(1) < W.at<float>(2) ||
+            W.at<float>(0) < W.at<float>(1) ||
+            norm(U*Mat::diag(W)*Vt, Q, CV_C) > FLT_EPSILON )
+            throw test_excep();
     }
     catch(const test_excep&)
     {
index 947d76d..6b9341b 100644 (file)
@@ -160,7 +160,7 @@ if(HAVE_FFMPEG)
     list(APPEND HIGHGUI_LIBRARIES ${BZIP2_LIBRARIES})
   endif()
   if(APPLE)
-    list(APPEND HIGHGUI_LIBRARIES "-framework VideoDecodeAcceleration")
+    list(APPEND HIGHGUI_LIBRARIES "-framework VideoDecodeAcceleration" bz2)
   endif()
 endif(HAVE_FFMPEG)