diff --git a/VectoCommon/VectoCommon/Utils/VectoMath.cs b/VectoCommon/VectoCommon/Utils/VectoMath.cs index 46c6642a64ab518ff54042c2346946d33fcf41e5..6d0684f7d7da9bd695cf4426331223d4f70fc6ee 100644 --- a/VectoCommon/VectoCommon/Utils/VectoMath.cs +++ b/VectoCommon/VectoCommon/Utils/VectoMath.cs @@ -378,6 +378,41 @@ namespace TUGraz.VectoCommon.Utils { return x < 0 ? -Math.Pow(-x, 1.0 / 3.0) : Math.Pow(x, 1.0 / 3.0); } + + + public static void LeastSquaresFitting<T>(IEnumerable<T> entries, Func<T, double> getX, Func<T, double> getY, + out double k, out double d, out double r) + { + // algoritm taken from http://mathworld.wolfram.com/LeastSquaresFitting.html (eqn. 27 & 28) + var count = 0; + var sumX = 0.0; + var sumY = 0.0; + var sumXSquare = 0.0; + var sumYSquare = 0.0; + var sumXY = 0.0; + foreach (var entry in entries) { + var x = getX(entry); + var y = getY(entry); + sumX += x; + sumY += y; + sumXSquare += x * x; + sumYSquare += y * y; + sumXY += x * y; + count++; + } + if (count == 0) { + k = 0; + d = 0; + r = 0; + return; + } + var ssxx = sumXSquare - sumX * sumX / count; + var ssxy = sumXY - sumX * sumY / count; + var ssyy = sumYSquare - sumY * sumY / count; + k = ssxy / ssxx; + d = (sumY - k * sumX) / count; + r = ssxy * ssxy / ssxx / ssyy; + } } [DebuggerDisplay("(X:{X}, Y:{Y}, Z:{Z})")] diff --git a/VectoCore/VectoCoreTest/Utils/VectoMathTest.cs b/VectoCore/VectoCoreTest/Utils/VectoMathTest.cs index 9c8c36324e545fb05f3844d21a494fa0eda4e90d..c9c082cb14365a0edb11158e4c85e0e9bb1156cb 100644 --- a/VectoCore/VectoCoreTest/Utils/VectoMathTest.cs +++ b/VectoCore/VectoCoreTest/Utils/VectoMathTest.cs @@ -104,5 +104,46 @@ namespace TUGraz.VectoCore.Tests.Utils Assert.AreEqual(0, cmp, 1e-15); } } + + [TestCase()] + public void TestLeastSquaresFittingExact() + { + var entries = new[] { + new { X = 0, Y = 4 }, + new { X = 10, Y = 8 }, + new { X = 20, Y = 12 } + }; + + double k, d, r; + VectoMath.LeastSquaresFitting(entries, x => x.X, x => x.Y, out k, out d, out r); + + Assert.AreEqual(4, d, 1e-6); + Assert.AreEqual(0.4, k, 1e-6); + Assert.AreEqual(1, r, 1e-6); + } + + [TestCase()] + public void TestLeastSquaresFittingEx1() + { + var entries = new[] { + new { X = 12, Y = 34.12 }, + new { X = 19, Y = 40.94 }, + new { X = 23, Y = 33.58 }, + new { X = 30, Y = 38.95 }, + new { X = 22, Y = 35.42}, + new { X = 11, Y = 32.12 }, + new { X = 13, Y = 28.57 }, + new { X = 28, Y = 40.97 }, + new { X = 10, Y = 32.06 }, + new { X = 11, Y = 30.55 }, + }; + + double k, d, r; + VectoMath.LeastSquaresFitting(entries, x => x.X, x => x.Y, out k, out d, out r); + + Assert.AreEqual(27.003529, d, 1e-6); + Assert.AreEqual(0.431535, k, 1e-6); + //Assert.AreEqual(1, r, 1e-3); + } } } \ No newline at end of file