From b75bac7975c525e7fda4dc3f861e99c8a250226d Mon Sep 17 00:00:00 2001
From: LaurentBerger <laurent.berger@univ-lemans.fr>
Date: Tue, 9 Aug 2016 16:00:03 +0200
Subject: [PATCH] Solve Issue 7063

consequences of changes

accuracy test

Solve issue 7063
---
 modules/core/src/matrix_decomp.cpp  |  8 ++++++-
 modules/core/test/test_math.cpp     | 19 ++++++++++++++++
 modules/ml/src/inner_functions.cpp  | 35 +++++------------------------
 modules/stitching/src/autocalib.cpp |  3 ---
 4 files changed, 32 insertions(+), 33 deletions(-)

diff --git a/modules/core/src/matrix_decomp.cpp b/modules/core/src/matrix_decomp.cpp
index dc9c014930..5b0d62a967 100644
--- a/modules/core/src/matrix_decomp.cpp
+++ b/modules/core/src/matrix_decomp.cpp
@@ -153,8 +153,12 @@ CholImpl(_Tp* A, size_t astep, int m, _Tp* b, size_t bstep, int n)
         L[i*astep + i] = (_Tp)(1./std::sqrt(s));
     }
 
-    if( !b )
+    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;
 }
diff --git a/modules/core/test/test_math.cpp b/modules/core/test/test_math.cpp
index 5eb9fe267f..32888ff734 100644
--- a/modules/core/test/test_math.cpp
+++ b/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. */
diff --git a/modules/ml/src/inner_functions.cpp b/modules/ml/src/inner_functions.cpp
index b39d403434..20abf7b052 100644
--- a/modules/ml/src/inner_functions.cpp
+++ b/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
diff --git a/modules/stitching/src/autocalib.cpp b/modules/stitching/src/autocalib.cpp
index 24145242c7..18b6e048d0 100644
--- a/modules/stitching/src/autocalib.cpp
+++ b/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;
 }