From 1db6a8a1f33786880706426313d72725f5f7bf99 Mon Sep 17 00:00:00 2001
From: Vincent Rabaud <vrabaud@google.com>
Date: Fri, 31 May 2024 19:51:58 +0200
Subject: [PATCH] Merge pull request #25665 from vrabaud:jacobian

Fix Homography computation. #25665

The bug was introduced in https://github.com/opencv/opencv/pull/25308

I am sorry I do not have a proper test.

### Pull Request Readiness Checklist

See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request

- [x] I agree to contribute to the project under Apache 2 License.
- [x] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV
- [x] The PR is proposed to the proper branch
- [x] There is a reference to the original bug report and related work
- [ ] There is accuracy test, performance test and test data in opencv_extra repository, if applicable
      Patch to opencv_extra has the same branch name.
- [ ] The feature is well documented and sample code can be built with the project CMake
---
 modules/calib3d/src/fundam.cpp           | 11 ++++----
 modules/calib3d/test/test_homography.cpp | 36 ++++++++++++++++++++++++
 2 files changed, 42 insertions(+), 5 deletions(-)

diff --git a/modules/calib3d/src/fundam.cpp b/modules/calib3d/src/fundam.cpp
index 599250114b..75d3456595 100644
--- a/modules/calib3d/src/fundam.cpp
+++ b/modules/calib3d/src/fundam.cpp
@@ -261,11 +261,11 @@ public:
             if( Jptr )
             {
                 Jptr[0] = Mx*ww; Jptr[1] = My*ww; Jptr[2] = ww;
-                Jptr[6] = -Mx*ww*xi; Jptr[7] = -My*ww*xi;
-                Jptr[11] = Mx*ww; Jptr[12] = My*ww; Jptr[13] = ww;
-                Jptr[14] = -Mx*ww*yi; Jptr[15] = -My*ww*yi;
+                Jptr[6] = -Mx*ww*xi; Jptr[7] = -My*ww*xi; Jptr[8] = -ww*xi;
+                Jptr[12] = Mx*ww; Jptr[13] = My*ww; Jptr[14] = ww;
+                Jptr[15] = -Mx*ww*yi; Jptr[16] = -My*ww*yi; Jptr[17] = -ww*yi;
 
-                Jptr += 16;
+                Jptr += 18;
             }
         }
 
@@ -274,7 +274,7 @@ public:
 
     Mat src, dst;
 };
-} // end namesapce cv
+} // end namespace cv
 
 namespace cv{
 static bool createAndRunRHORegistrator(double confidence,
@@ -426,6 +426,7 @@ cv::Mat cv::findHomography( InputArray _points1, InputArray _points2,
                 cb->runKernel( src, dst, H );
             Mat H8(9, 1, CV_64F, H.ptr<double>());
             LMSolver::create(makePtr<HomographyRefineCallback>(src, dst), 10)->run(H8);
+            H.convertTo(H, H.type(), scaleFor(H.at<double>(2,2)));
         }
     }
 
diff --git a/modules/calib3d/test/test_homography.cpp b/modules/calib3d/test/test_homography.cpp
index 6a9cdcdd72..d8ab243a84 100644
--- a/modules/calib3d/test/test_homography.cpp
+++ b/modules/calib3d/test/test_homography.cpp
@@ -727,4 +727,40 @@ TEST(Calib3d_Homography, not_normalized)
     }
 }
 
+TEST(Calib3d_Homography, Refine)
+{
+    Mat_<double> p1({10, 2}, {41, -86, -87, 99, 66, -96, -86, -8, -67, 24,
+                              -87, -76, -19, 89, 37, -4, -86, -86, -66, -53});
+    Mat_<double> p2({10, 2}, {
+        0.007723226608700208, -1.177541410622515,
+        -0.1909072353027552, -0.4247610181930323,
+        -0.134992319993638, -0.6469949816560389,
+        -0.3570627451405215, 0.1811469436293486,
+        -0.3005671881038939, -0.02325733734262935,
+        -0.4404509481789249, 0.4851526464158342,
+        0.6343346428859541, -3.396187657072353,
+        -0.3539383967092603, 0.1469447227353143,
+        -0.4526924606856586, 0.5296757109061794,
+        -0.4309974583614644, 0.4522732662733471
+    });
+    hconcat(p1, Mat::ones(p1.rows, 1, CV_64F), p1);
+    hconcat(p2, Mat::ones(p2.rows, 1, CV_64F), p2);
+
+    for(int method : std::vector<int>({0, RANSAC, LMEDS}))
+    {
+        Mat h = findHomography(p1, p2, method);
+        EXPECT_NEAR(h.at<double>(2, 2), 1.0, 1e-7);
+
+        Mat proj = p1 * h.t();
+        proj.col(0) /= proj.col(2);
+        proj.col(1) /= proj.col(2);
+
+        Mat error;
+        cv::pow(p2.colRange(0, 2) - proj.colRange(0, 2), 2, error);
+        cv::reduce(error, error, 1, REDUCE_SUM);
+        cv::reduce(error, error, 0, REDUCE_AVG);
+        EXPECT_LE(sqrt(error.at<double>(0, 0)), method == LMEDS ? 7e-4 : 7e-5);
+    }
+}
+
 }} // namespace