Merge pull request #7077 from LaurentBerger:I7063

modules/core/src/matrix_decomp.cpp

@@ -154,7 +154,11 @@ CholImpl(_Tp* A, size_t astep, int m, _Tp* b, size_t bstep, int n)
     }
 
     if (!b)
+    {
+        for( i = 0; i < m; i++ )
+             L[i*astep + i]=1/L[i*astep + i];
         return true;
+   }
 
     // LLt x = b
     // 1: L y = b
@@ -193,6 +197,8 @@ CholImpl(_Tp* A, size_t astep, int m, _Tp* b, size_t bstep, int n)
             b[i*bstep + j] = (_Tp)(s*L[i*astep + i]);
         }
     }
+    for( i = 0; i < m; i++ )
+            L[i*astep + i]=1/L[i*astep + i];
 
     return true;
 }

modules/core/test/test_math.cpp

@@ -2977,4 +2977,23 @@ TEST(Core_Pow, special)
     }
 }
 
+TEST(Core_Cholesky, accuracy64f)
+{
+    const int n = 5;
+    Mat A(n, n, CV_64F), refA;
+    Mat mean(1, 1, CV_64F);
+    *mean.ptr<double>() = 10.0;
+    Mat dev(1, 1, CV_64F);
+    *dev.ptr<double>() = 10.0;
+    RNG rng(10);
+    rng.fill(A, RNG::NORMAL, mean, dev);
+    A = A*A.t();
+    A.copyTo(refA);
+    Cholesky(A.ptr<double>(), A.step, n, NULL, 0, 0);
+
+   for (int i = 0; i < A.rows; i++)
+       for (int j = i + 1; j < A.cols; j++)
+           A.at<double>(i, j) = 0.0;
+   EXPECT_TRUE(norm(refA - A*A.t()) < 10e-5);
+}
 /* End of file. */

modules/ml/src/inner_functions.cpp

@@ -116,35 +116,12 @@ static void Cholesky( const Mat& A, Mat& S )
 {
     CV_Assert(A.type() == CV_32F);
 
-    int dim = A.rows;
-    S.create(dim, dim, CV_32F);
-
-    int i, j, k;
-
-    for( i = 0; i < dim; i++ )
-    {
-        for( j = 0; j < i; j++ )
-            S.at<float>(i,j) = 0.f;
-
-        float sum = 0.f;
-        for( k = 0; k < i; k++ )
-        {
-            float val = S.at<float>(k,i);
-            sum += val*val;
-        }
-
-        S.at<float>(i,i) = std::sqrt(std::max(A.at<float>(i,i) - sum, 0.f));
-        float ival = 1.f/S.at<float>(i, i);
-
-        for( j = i + 1; j < dim; j++ )
-        {
-            sum = 0;
-            for( k = 0; k < i; k++ )
-                sum += S.at<float>(k, i) * S.at<float>(k, j);
-
-            S.at<float>(i, j) = (A.at<float>(i, j) - sum)*ival;
-        }
-    }
+    S = A.clone();
+    cv::Cholesky ((float*)S.ptr(),S.step, S.rows,NULL, 0, 0);
+    S = S.t();
+    for (int i=1;i<S.rows;i++)
+        for (int j=0;j<i;j++)
+            S.at<float>(i,j)=0;
 }
 
 /* Generates <sample> from multivariate normal distribution, where <mean> - is an

modules/stitching/src/autocalib.cpp

@@ -51,9 +51,6 @@ static inline bool decomposeCholesky(double* A, size_t astep, int m)
 {
     if (!hal::Cholesky64f(A, astep, m, 0, 0, 0))
         return false;
-    astep /= sizeof(A[0]);
-    for (int i = 0; i < m; ++i)
-        A[i*astep + i] = (double)(1./A[i*astep + i]);
     return true;
 }