improved accuracy of the matrix determinant and matrix inversion functions (trac #431)

modules/core/src/lapack.cpp

@@ -229,11 +229,13 @@ double determinant( const Mat& mat )
     size_t step = mat.step;
     const uchar* m = mat.data;
 
-    CV_Assert( mat.rows == mat.cols );
+    CV_Assert( mat.rows == mat.cols && (type == CV_32F || type == CV_64F));
 
     #define Mf(y, x) ((float*)(m + y*step))[x]
     #define Md(y, x) ((double*)(m + y*step))[x]
 
+    if( rows <= 10 )
+    {
         if( type == CV_32F )
         {
             if( rows == 2 )
@@ -244,30 +246,21 @@ double determinant( const Mat& mat )
                 result = Mf(0,0);
             else
             {
-            integer i, n = rows, *ipiv, info=0;
-            int bufSize = n*n*sizeof(float) + (n+1)*sizeof(ipiv[0]), sign=0;
+                size_t bufSize = rows*rows*sizeof(float);
                 AutoBuffer<uchar> buffer(bufSize);
-
-            Mat a(n, n, CV_32F, (uchar*)buffer);
+                Mat a(rows, rows, CV_32F, (uchar*)buffer);
                 mat.copyTo(a);
                 
-            ipiv = (integer*)cvAlignPtr(a.data + a.step*a.rows, sizeof(integer));
-            sgetrf_(&n, &n, (float*)a.data, &n, ipiv, &info);
-            assert(info >= 0);
-
-            if( info == 0 )
+                result = LU((float*)a.data, rows, 0, 0);
+                if( result )
                 {
-                result = 1;
-                for( i = 0; i < n; i++ )
-                {
-                    result *= ((float*)a.data)[i*(n+1)];
-                    sign ^= ipiv[i] != i+1;
-                }
-                result *= sign ? -1 : 1;
+                    for( int i = 0; i < rows; i++ )
+                        result *= ((const float*)(a.data + a.step*i))[i];
+                    result = 1./result;
                 }
             }
         }
-    else if( type == CV_64F )
+        else
         {
             if( rows == 2 )
                 result = det2(Md);
@@ -277,14 +270,31 @@ double determinant( const Mat& mat )
                 result = Md(0,0);
             else
             {
-            integer i, n = rows, *ipiv, info=0;
-            int bufSize = n*n*sizeof(double) + (n+1)*sizeof(ipiv[0]), sign=0;
+                size_t bufSize = rows*rows*sizeof(double);
+                AutoBuffer<uchar> buffer(bufSize);
+                Mat a(rows, rows, CV_64F, (uchar*)buffer);
+                mat.copyTo(a);
+                
+                result = LU((double*)a.data, rows, 0, 0);
+                if( result )
+                {
+                    for( int i = 0; i < rows; i++ )
+                        result *= ((const double*)(a.data + a.step*i))[i];
+                    result = 1./result;
+                }
+            }
+        }
+    }
+    else
+    {
+        integer i, n = rows, *ipiv, info=0, sign = 0;
+        size_t bufSize = n*n*sizeof(double) + (n+1)*sizeof(ipiv[0]);
         AutoBuffer<uchar> buffer(bufSize);
 
         Mat a(n, n, CV_64F, (uchar*)buffer);
-            mat.copyTo(a);
-            ipiv = (integer*)cvAlignPtr(a.data + a.step*a.rows, sizeof(integer));
+        mat.convertTo(a, CV_64F);
 
+        ipiv = (integer*)cvAlignPtr(a.data + a.step*a.rows, sizeof(integer));
         dgetrf_(&n, &n, (double*)a.data, &n, ipiv, &info);
         assert(info >= 0);
 
@@ -299,9 +309,6 @@ double determinant( const Mat& mat )
             result *= sign ? -1 : 1;
         }
     }
-    }
-    else
-        CV_Error( CV_StsUnsupportedFormat, "" );
 
     #undef Mf
     #undef Md
@@ -341,8 +348,6 @@ double invert( const Mat& src, Mat& dst, int method )
     CV_Assert( src.rows == src.cols && (type == CV_32F || type == CV_64F));
     dst.create( src.rows, src.cols, type );
 
-    if( method == DECOMP_LU || method == DECOMP_CHOLESKY )
-    {
     if( src.rows <= 3 )
     {
         uchar* srcdata = src.data;
@@ -470,64 +475,89 @@ double invert( const Mat& src, Mat& dst, int method )
         return result;
     }
 
-        src.copyTo(dst);
-        integer n = dst.cols, lwork=-1, elem_size = CV_ELEM_SIZE(type),
-            lda = (int)(dst.step/elem_size), piv1=0, info=0;    
-            
-        if( method == DECOMP_LU )
+    if( dst.cols <= 10 )
     {
-            int buf_size = (int)(n*sizeof(integer));
-            AutoBuffer<uchar> buf;
-            uchar* buffer;
+        int n = dst.cols, elem_size = CV_ELEM_SIZE(type);
+        AutoBuffer<uchar> buf(n*n*2*elem_size);
+        Mat src1(n, n, type, (uchar*)buf);
+        Mat dst1(n, n, type, dst.isContinuous() ? dst.data : src1.data + n*n*elem_size);
+        src.copyTo(src1);
+        setIdentity(dst1);
         
-            if( type == CV_32F )
-            {
-                real work1 = 0;
-                sgetri_(&n, (float*)dst.data, &lda, &piv1, &work1, &lwork, &info);
-                lwork = cvRound(work1);
+        if( method == DECOMP_LU && type == CV_32F )
+            result = LU((float*)src1.data, n, (float*)dst1.data, n);
+        else if( method == DECOMP_LU && type == CV_64F )
+            result = LU((double*)src1.data, n, (double*)dst1.data, n);
+        else if( method == DECOMP_LU && type == CV_32F )
+            result = Cholesky((float*)src1.data, n, (float*)dst1.data, n);
+        else
+            result = Cholesky((double*)src1.data, n, (double*)dst1.data, n);
+        dst1.copyTo(dst);
+        result = std::abs(result);
     }
     else
+    {
+        integer n = dst.cols, lwork=-1, lda = n, piv1=0, info=0;
+        int t_size = type == CV_32F ? n*n*sizeof(double) : 0;
+        int buf_size = t_size;
+        AutoBuffer<uchar> buf;
+        
+        if( method == DECOMP_LU )
         {
             double work1 = 0;
             dgetri_(&n, (double*)dst.data, &lda, &piv1, &work1, &lwork, &info);
             lwork = cvRound(work1);
-            }
 
-            buf_size += (int)((lwork + 1)*elem_size);
+            buf_size += (int)(n*sizeof(integer) + (lwork + 1)*sizeof(double));
             buf.allocate(buf_size);
-            buffer = (uchar*)buf;
+            uchar* buffer = (uchar*)buf;
             
+            Mat arr = dst;
             if( type == CV_32F )
             {
-                sgetrf_(&n, &n, (float*)dst.data, &lda, (integer*)buffer, &info);
-                if(info==0)
-                    sgetri_(&n, (float*)dst.data, &lda, (integer*)buffer,
-                        (float*)(buffer + n*sizeof(integer)), &lwork, &info);
+                arr = Mat(n, n, CV_64F, buffer);
+                src.convertTo(arr, CV_64F);
+                buffer += t_size;
             }
             else
             {
-                dgetrf_(&n, &n, (double*)dst.data, &lda, (integer*)buffer, &info);
-                if(info==0)
-                    dgetri_(&n, (double*)dst.data, &lda, (integer*)buffer,
-                        (double*)cvAlignPtr(buffer + n*sizeof(integer), elem_size), &lwork, &info);
+                src.copyTo(arr);
+                lda = arr.step/sizeof(double);
             }
+
+            dgetrf_(&n, &n, (double*)arr.data, &lda, (integer*)buffer, &info);
+            if(info==0)
+                dgetri_(&n, (double*)arr.data, &lda, (integer*)buffer,
+                    (double*)cvAlignPtr(buffer + n*sizeof(integer), sizeof(double)),
+                    &lwork, &info);
+            if(info==0 && arr.data != dst.data)
+                arr.convertTo(dst, dst.type());
         }
-        else if( method == CV_CHOLESKY )
+        else if( method == DECOMP_CHOLESKY )
         {
-            char L[] = {'L', '\0'};
+            Mat arr = dst;
             if( type == CV_32F )
             {
-                spotrf_(L, &n, (float*)dst.data, &lda, &info);
-                if(info==0)
-                    spotri_(L, &n, (float*)dst.data, &lda, &info);
+                buf.allocate(buf_size);
+                arr = Mat(n, n, CV_64F, (uchar*)buf);
+                src.convertTo(arr, CV_64F);
             }
             else
             {
-                dpotrf_(L, &n, (double*)dst.data, &lda, &info);
+                src.copyTo(arr);
+                lda = arr.step/sizeof(double);
+            }
+            
+            char L[] = {'L', '\0'};
+            dpotrf_(L, &n, (double*)arr.data, &lda, &info);
+            if(info==0)
+                dpotri_(L, &n, (double*)arr.data, &lda, &info);
             if(info==0)
-                    dpotri_(L, &n, (double*)dst.data, &lda, &info);
+            {
+                completeSymm(arr);
+                if( arr.data != dst.data )
+                    arr.convertTo(dst, dst.type());
             }
-            completeSymm(dst);
         }
         result = info == 0;
     }

tests/cxcore/src/amath.cpp

@@ -2420,7 +2420,7 @@ void CxCore_DetTest::prepare_to_validation( int )
     *((CvScalar*)(test_mat[REF_OUTPUT][0].data.db)) = cvRealScalar(cvTsLU(&test_mat[TEMP][0], 0, 0));
 }
 
-//CxCore_DetTest det_test;
+CxCore_DetTest det_test;
 
 
 
@@ -2475,8 +2475,8 @@ void CxCore_InvertTest::get_test_array_types_and_sizes( int test_case_idx, CvSiz
         if( bits & 4 )
         {
             sizes[INPUT][0] = cvSize(min_size, min_size);
-            if( bits & 8 )
-                method = CV_SVD_SYM;
+            if( bits & 16 )
+                method = CV_CHOLESKY;
         }
     }
     else
@@ -2536,7 +2536,7 @@ int CxCore_InvertTest::prepare_test_case( int test_case_idx )
     {
         cvTsFloodWithZeros( &test_mat[INPUT][0], ts->get_rng() );
 
-        if( method == CV_SVD_SYM )
+        if( method == CV_CHOLESKY )
         {
             cvTsGEMM( &test_mat[INPUT][0], &test_mat[INPUT][0], 1.,
                       0, 0., &test_mat[TEMP][0], CV_GEMM_B_T );