From 0e019e3a43bd3859ac51f8a66d6f7d63c153a139 Mon Sep 17 00:00:00 2001
From: wmayer <wmayer@users.sourceforge.net>
Date: Wed, 6 May 2020 23:25:18 +0200
Subject: [PATCH] Mesh: [skip ci] improve cylinder fit

---
 src/Mod/Mesh/App/Core/Approximation.cpp | 31 +++++++++++++++++++++++--
 src/Mod/Mesh/App/Core/CylinderFit.cpp   | 16 ++++++++-----
 2 files changed, 39 insertions(+), 8 deletions(-)

diff --git a/src/Mod/Mesh/App/Core/Approximation.cpp b/src/Mod/Mesh/App/Core/Approximation.cpp
index f811f32f62..ee7877abe9 100644
--- a/src/Mod/Mesh/App/Core/Approximation.cpp
+++ b/src/Mod/Mesh/App/Core/Approximation.cpp
@@ -1037,6 +1037,32 @@ CylinderFit::~CylinderFit()
 
 Base::Vector3f CylinderFit::GetInitialAxisFromNormals(const std::vector<Base::Vector3f>& n) const
 {
+    int nc = 0;
+    double x = 0.0;
+    double y = 0.0;
+    double z = 0.0;
+    for (int i = 0; i < (int)n.size()-1; ++i) {
+        for (int j = i+1; j < (int)n.size(); ++j) {
+            Base::Vector3f cross = n[i] % n[j];
+            if (cross.Sqr() > 1.0e-6) {
+                cross.Normalize();
+                x += cross.x;
+                y += cross.y;
+                z += cross.z;
+                ++nc;
+            }
+        }
+    }
+
+    if (nc > 0) {
+        x /= (double)nc;
+        y /= (double)nc;
+        z /= (double)nc;
+        Base::Vector3f axis(x,y,z);
+        axis.Normalize();
+        return axis;
+    }
+
     PlaneFit planeFit;
     planeFit.AddPoints(n);
     planeFit.Fit();
@@ -1082,11 +1108,12 @@ float CylinderFit::Fit()
         _vBase.x, _vBase.y, _vBase.z, _vAxis.x, _vAxis.y, _vAxis.z, _fRadius, GetStdDeviation());
 #endif
 
+    // Do the cylinder fit
     MeshCoreFit::CylinderFit cylFit;
     cylFit.AddPoints(_vPoints);
-    //cylFit.SetApproximations(_fRadius, Base::Vector3d(_vBase.x, _vBase.y, _vBase.z), Base::Vector3d(_vAxis.x, _vAxis.y, _vAxis.z));
+    if (_fLastResult < FLOAT_MAX)
+        cylFit.SetApproximations(_fRadius, Base::Vector3d(_vBase.x, _vBase.y, _vBase.z), Base::Vector3d(_vAxis.x, _vAxis.y, _vAxis.z));
 
-    // Do the cylinder fit
     float result = cylFit.Fit();
     if (result < FLOAT_MAX) {
         Base::Vector3d base = cylFit.GetBase();
diff --git a/src/Mod/Mesh/App/Core/CylinderFit.cpp b/src/Mod/Mesh/App/Core/CylinderFit.cpp
index b66e0470bd..d5bfb3480f 100644
--- a/src/Mod/Mesh/App/Core/CylinderFit.cpp
+++ b/src/Mod/Mesh/App/Core/CylinderFit.cpp
@@ -359,26 +359,30 @@ void CylinderFit::findBestSolDirection(SolutionD &solDir)
 	if (biggest < 0.0)
 		dir.Set(-dir.x, -dir.y, -dir.z);	// multiplies by -1
 
+	double fixedVal = 0.0;
 	double lambda;
 	switch (solDir)
 	{
 		case solL:
-			lambda = -pos.x / dir.x;
-			pos.x = 0.0;//meanXObs();
+			fixedVal = meanXObs();
+			lambda = (fixedVal - pos.x) / dir.x;
+			pos.x = fixedVal;
 			pos.y = pos.y + lambda * dir.y;
 			pos.z = pos.z + lambda * dir.z;
 			break;
 		case solM:
-			lambda = -pos.y / dir.y;
+			fixedVal = meanYObs();
+			lambda = (fixedVal - pos.y) / dir.y;
 			pos.x = pos.x + lambda * dir.x;
-			pos.y = 0.0;//meanYObs();
+			pos.y = fixedVal;
 			pos.z = pos.z + lambda * dir.z;
 			break;
 		case solN:
-			lambda = -pos.z / dir.z;
+			fixedVal = meanZObs();
+			lambda = (fixedVal - pos.z) / dir.z;
 			pos.x = pos.x + lambda * dir.x;
 			pos.y = pos.y + lambda * dir.y;
-			pos.z = 0.0;//meanZObs();
+			pos.z = fixedVal;
 			break;
 	}
 	_vAxis = dir;