Use double precision for accuracy

src/ImageSharp/Common/Helpers/GaussianEliminationSolver.cs

@@ -20,24 +20,24 @@ internal static class GaussianEliminationSolver
     /// The matrix passed to this method must be a square matrix.
     /// If the matrix is singular (i.e., has no unique solution), an <see cref="NotSupportedException"/> will be thrown.
     /// </remarks>
-    public static void Solve(float[][] matrix, float[] result)
+    public static void Solve(double[][] matrix, double[] result)
     {
         TransformToRowEchelonForm(matrix, result);
         ApplyBackSubstitution(matrix, result);
     }
 
-    private static void TransformToRowEchelonForm(float[][] matrix, float[] result)
+    private static void TransformToRowEchelonForm(double[][] matrix, double[] result)
     {
         int colCount = matrix.Length;
         int rowCount = matrix[0].Length;
         int pivotRow = 0;
         for (int pivotCol = 0; pivotCol < colCount; pivotCol++)
         {
-            float maxValue = float.Abs(matrix[pivotRow][pivotCol]);
+            double maxValue = double.Abs(matrix[pivotRow][pivotCol]);
             int maxIndex = pivotRow;
             for (int r = pivotRow + 1; r < rowCount; r++)
             {
-                float value = float.Abs(matrix[r][pivotCol]);
+                double value = double.Abs(matrix[r][pivotCol]);
                 if (value > maxValue)
                 {
                     maxIndex = r;
@@ -55,7 +55,7 @@ internal static class GaussianEliminationSolver
 
             for (int r = pivotRow + 1; r < rowCount; r++)
             {
-                float fraction = matrix[r][pivotCol] / matrix[pivotRow][pivotCol];
+                double fraction = matrix[r][pivotCol] / matrix[pivotRow][pivotCol];
                 for (int c = pivotCol + 1; c < colCount; c++)
                 {
                     matrix[r][c] -= matrix[pivotRow][c] * fraction;
@@ -69,7 +69,7 @@ internal static class GaussianEliminationSolver
         }
     }
 
-    private static void ApplyBackSubstitution(float[][] matrix, float[] result)
+    private static void ApplyBackSubstitution(double[][] matrix, double[] result)
     {
         int rowCount = matrix[0].Length;
 

src/ImageSharp/Common/Helpers/QuadDistortionHelper.cs

@@ -40,7 +40,7 @@ internal static class QuadDistortionHelper
         PointF q3 = bottomRight;
         PointF q4 = bottomLeft;
 
-        float[][] matrixData =
+        double[][] matrixData =
         [
             [p1.X, p1.Y, 1, 0, 0, 0, -p1.X * q1.X, -p1.Y * q1.X],
             [0, 0, 0, p1.X, p1.Y, 1, -p1.X * q1.Y, -p1.Y * q1.Y],
@@ -52,7 +52,7 @@ internal static class QuadDistortionHelper
             [0, 0, 0, p4.X, p4.Y, 1, -p4.X * q4.Y, -p4.Y * q4.Y],
         ];
 
-        float[] b =
+        double[] b =
         [
             q1.X,
             q1.Y,
@@ -68,10 +68,10 @@ internal static class QuadDistortionHelper
 
 #pragma warning disable SA1117
         Matrix4x4 projectionMatrix = new(
-            b[0], b[3], 0, b[6],
-            b[1], b[4], 0, b[7],
+            (float)b[0], (float)b[3], 0, (float)b[6],
+            (float)b[1], (float)b[4], 0, (float)b[7],
             0, 0, 1, 0,
-            b[2], b[5], 0, 1);
+            (float)b[2], (float)b[5], 0, 1);
 #pragma warning restore SA1117
 
         return projectionMatrix;

tests/ImageSharp.Tests/Common/GaussianEliminationSolverTest.cs

@@ -9,7 +9,7 @@ public class GaussianEliminationSolverTest
 {
     [Theory]
     [MemberData(nameof(MatrixTestData))]
-    public void CanSolve(float[][] matrix, float[] result, float[] expected)
+    public void CanSolve(double[][] matrix, double[] result, double[] expected)
     {
         GaussianEliminationSolver.Solve(matrix, result);
 
@@ -19,44 +19,44 @@ public class GaussianEliminationSolverTest
         }
     }
 
-    public static TheoryData<float[][], float[], float[]> MatrixTestData
+    public static TheoryData<double[][], double[], double[]> MatrixTestData
     {
         get
         {
-            TheoryData<float[][], float[], float[]> data = [];
+            TheoryData<double[][], double[], double[]> data = [];
             {
-                float[][] matrix =
+                double[][] matrix =
                 [
                     [2, 3, 4],
                     [1, 2, 3],
                     [3, -4, 0],
                 ];
-                float[] result = [6, 4, 10];
-                float[] expected = [18 / 11f, -14 / 11f, 18 / 11f];
+                double[] result = [6, 4, 10];
+                double[] expected = [18 / 11f, -14 / 11f, 18 / 11f];
                 data.Add(matrix, result, expected);
             }
 
             {
-                float[][] matrix =
+                double[][] matrix =
                 [
                     [1, 4, -1],
                     [2, 5, 8],
                     [1, 3, -3],
                 ];
-                float[] result = [4, 15, 1];
-                float[] expected = [1, 1, 1];
+                double[] result = [4, 15, 1];
+                double[] expected = [1, 1, 1];
                 data.Add(matrix, result, expected);
             }
 
             {
-                float[][] matrix =
+                double[][] matrix =
                 [
                     [-1, 0, 0],
                     [0, 1, 0],
                     [0, 0, 1],
                 ];
-                float[] result = [1, 2, 3];
-                float[] expected = [-1, 2, 3];
+                double[] result = [1, 2, 3];
+                double[] expected = [-1, 2, 3];
                 data.Add(matrix, result, expected);
             }