Add QuadDistortion to ProjectiveTransformBuilder

Author: Socolin

src/ImageSharp/Common/Helpers/GaussianEliminationSolver.cs

@@ -0,0 +1,89 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+using System.Numerics;
+
+namespace SixLabors.ImageSharp.Common.Helpers;
+
+/// <summary>
+/// Represents a solver for systems of linear equations using the Gaussian Elimination method.
+/// This class applies Gaussian Elimination to transform the matrix into row echelon form and then performs back substitution to find the solution vector.
+/// This implementation is based on: https://www.algorithm-archive.org/contents/gaussian_elimination/gaussian_elimination.html
+/// </summary>
+/// <typeparam name="TNumber">The type of numbers used in the matrix and solution vector. Must implement the <see cref="INumber{TNumber}"/> interface.</typeparam>
+internal static class GaussianEliminationSolver<TNumber>
+    where TNumber : INumber<TNumber>
+{
+    /// <summary>
+    /// Solves the system of linear equations represented by the given matrix and result vector using Gaussian Elimination.
+    /// </summary>
+    /// <param name="matrix">The square matrix representing the coefficients of the linear equations.</param>
+    /// <param name="result">The vector representing the constants on the right-hand side of the linear equations.</param>
+    /// <exception cref="Exception">Thrown if the matrix is singular and cannot be solved.</exception>
+    /// <remarks>
+    /// 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(TNumber[][] matrix, TNumber[] result)
+    {
+        TransformToRowEchelonForm(matrix, result);
+        ApplyBackSubstitution(matrix, result);
+    }
+
+    private static void TransformToRowEchelonForm(TNumber[][] matrix, TNumber[] result)
+    {
+        int colCount = matrix.Length;
+        int rowCount = matrix[0].Length;
+        int pivotRow = 0;
+        for (int pivotCol = 0; pivotCol < colCount; pivotCol++)
+        {
+            TNumber maxValue = TNumber.Abs(matrix[pivotRow][pivotCol]);
+            int maxIndex = pivotRow;
+            for (int r = pivotRow + 1; r < rowCount; r++)
+            {
+                TNumber value = TNumber.Abs(matrix[r][pivotCol]);
+                if (value > maxValue)
+                {
+                    maxIndex = r;
+                    maxValue = value;
+                }
+            }
+
+            if (matrix[maxIndex][pivotCol] == TNumber.Zero)
+            {
+                throw new NotSupportedException("Matrix is singular and cannot be solve");
+            }
+
+            (matrix[pivotRow], matrix[maxIndex]) = (matrix[maxIndex], matrix[pivotRow]);
+            (result[pivotRow], result[maxIndex]) = (result[maxIndex], result[pivotRow]);
+
+            for (int r = pivotRow + 1; r < rowCount; r++)
+            {
+                TNumber fraction = matrix[r][pivotCol] / matrix[pivotRow][pivotCol];
+                for (int c = pivotCol + 1; c < colCount; c++)
+                {
+                    matrix[r][c] -= matrix[pivotRow][c] * fraction;
+                }
+
+                result[r] -= result[pivotRow] * fraction;
+                matrix[r][pivotCol] = TNumber.Zero;
+            }
+
+            pivotRow++;
+        }
+    }
+
+    private static void ApplyBackSubstitution(TNumber[][] matrix, TNumber[] result)
+    {
+        int rowCount = matrix[0].Length;
+
+        for (int row = rowCount - 1; row >= 0; row--)
+        {
+            result[row] /= matrix[row][row];
+
+            for (int r = 0; r < row; r++)
+            {
+                result[r] -= result[row] * matrix[r][row];
+            }
+        }
+    }
+}

src/ImageSharp/Common/Helpers/QuadDistortionHelper.cs

@@ -0,0 +1,76 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using System.Numerics;
+
+namespace SixLabors.ImageSharp.Common.Helpers;
+
+/// <summary>
+/// Provides helper methods for performing quad distortion transformations.
+/// </summary>
+internal static class QuadDistortionHelper
+{
+    /// <summary>
+    /// Computes the projection matrix for a quad distortion transformation.
+    /// </summary>
+    /// <param name="rectangle">The source rectangle.</param>
+    /// <param name="topLeft">The top-left point of the distorted quad.</param>
+    /// <param name="topRight">The top-right point of the distorted quad.</param>
+    /// <param name="bottomRight">The bottom-right point of the distorted quad.</param>
+    /// <param name="bottomLeft">The bottom-left point of the distorted quad.</param>
+    /// <returns>The computed projection matrix for the quad distortion.</returns>
+    /// <remarks>
+    /// This method is based on the algorithm described in the following article:
+    /// https://blog.mbedded.ninja/mathematics/geometry/projective-transformations/
+    /// </remarks>
+    public static Matrix4x4 ComputeQuadDistortMatrix(Rectangle rectangle, PointF topLeft, PointF topRight, PointF bottomRight, PointF bottomLeft)
+    {
+        PointF p1 = new(rectangle.X, rectangle.Y);
+        PointF p2 = new(rectangle.X + rectangle.Width, rectangle.Y);
+        PointF p3 = new(rectangle.X + rectangle.Width, rectangle.Y + rectangle.Height);
+        PointF p4 = new(rectangle.X, rectangle.Y + rectangle.Height);
+
+        PointF q1 = topLeft;
+        PointF q2 = topRight;
+        PointF q3 = bottomRight;
+        PointF q4 = bottomLeft;
+
+        // @formatter:off
+        float[][] 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],
+            [p2.X, p2.Y, 1, 0,    0,    0, -p2.X * q2.X, -p2.Y * q2.X],
+            [0,    0,    0, p2.X, p2.Y, 1, -p2.X * q2.Y, -p2.Y * q2.Y],
+            [p3.X, p3.Y, 1, 0,    0,    0, -p3.X * q3.X, -p3.Y * q3.X],
+            [0,    0,    0, p3.X, p3.Y, 1, -p3.X * q3.Y, -p3.Y * q3.Y],
+            [p4.X, p4.Y, 1, 0,    0,    0, -p4.X * q4.X, -p4.Y * q4.X],
+            [0,    0,    0, p4.X, p4.Y, 1, -p4.X * q4.Y, -p4.Y * q4.Y],
+        ];
+
+        float[] b =
+        [
+            q1.X,
+            q1.Y,
+            q2.X,
+            q2.Y,
+            q3.X,
+            q3.Y,
+            q4.X,
+            q4.Y,
+        ];
+
+        GaussianEliminationSolver<float>.Solve(matrixData, b);
+
+ #pragma warning disable SA1117
+        Matrix4x4 projectionMatrix = new(
+            b[0], b[3], 0, b[6],
+            b[1], b[4], 0, b[7],
+            0,    0,    1, 0,
+            b[2], b[5], 0, 1);
+ #pragma warning restore SA1117
+
+        // @formatter:on
+        return projectionMatrix;
+    }
+}

src/ImageSharp/Processing/ProjectiveTransformBuilder.cs

@@ -2,6 +2,7 @@
 // Licensed under the Six Labors Split License.
 
 using System.Numerics;
+using SixLabors.ImageSharp.Common.Helpers;
 using SixLabors.ImageSharp.Processing.Processors.Transforms;
 
 namespace SixLabors.ImageSharp.Processing;
@@ -270,6 +271,28 @@ public ProjectiveTransformBuilder AppendTranslation(PointF position)
     public ProjectiveTransformBuilder AppendTranslation(Vector2 position)
         => this.AppendMatrix(Matrix4x4.CreateTranslation(new Vector3(position, 0)));
 
+    /// <summary>
+    /// Prepends a quad distortion matrix using the specified corner points.
+    /// </summary>
+    /// <param name="topLeft">The top-left corner point of the distorted quad.</param>
+    /// <param name="topRight">The top-right corner point of the distorted quad.</param>
+    /// <param name="bottomRight">The bottom-right corner point of the distorted quad.</param>
+    /// <param name="bottomLeft">The bottom-left corner point of the distorted quad.</param>
+    /// <returns>The <see cref="ProjectiveTransformBuilder"/>.</returns>
+    public ProjectiveTransformBuilder PrependQuadDistortion(PointF topLeft, PointF topRight, PointF bottomRight, PointF bottomLeft)
+        => this.Prepend(size => QuadDistortionHelper.ComputeQuadDistortMatrix(new Rectangle(Point.Empty, size), topLeft, topRight, bottomRight, bottomLeft));
+
+    /// <summary>
+    /// Appends a quad distortion matrix using the specified corner points.
+    /// </summary>
+    /// <param name="topLeft">The top-left corner point of the distorted quad.</param>
+    /// <param name="topRight">The top-right corner point of the distorted quad.</param>
+    /// <param name="bottomRight">The bottom-right corner point of the distorted quad.</param>
+    /// <param name="bottomLeft">The bottom-left corner point of the distorted quad.</param>
+    /// <returns>The <see cref="ProjectiveTransformBuilder"/>.</returns>
+    public ProjectiveTransformBuilder AppendQuadDistortion(PointF topLeft, PointF topRight, PointF bottomRight, PointF bottomLeft)
+        => this.Append(size => QuadDistortionHelper.ComputeQuadDistortMatrix(new Rectangle(Point.Empty, size), topLeft, topRight, bottomRight, bottomLeft));
+
     /// <summary>
     /// Prepends a raw matrix.
     /// </summary>

tests/ImageSharp.Tests/Common/GaussianEliminationSolverTest.cs

@@ -0,0 +1,66 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Common.Helpers;
+
+namespace SixLabors.ImageSharp.Tests.Common;
+
+public class GaussianEliminationSolverTest
+{
+    [Theory]
+    [MemberData(nameof(MatrixTestData))]
+    public void CanSolve(float[][] matrix, float[] result, float[] expected)
+    {
+        GaussianEliminationSolver<float>.Solve(matrix, result);
+
+        for (int i = 0; i < expected.Length; i++)
+        {
+            Assert.Equal(result[i], expected[i], 4);
+        }
+    }
+
+    public static TheoryData<float[][], float[], float[]> MatrixTestData
+    {
+        get
+        {
+            TheoryData<float[][], float[], float[]> data = [];
+            {
+                float[][] matrix =
+                [
+                    [2, 3, 4],
+                    [1, 2, 3],
+                    [3, -4, 0],
+                ];
+                float[] result = [6, 4, 10];
+                float[] expected = [18 / 11f, -14 / 11f, 18 / 11f];
+                data.Add(matrix, result, expected);
+            }
+
+            {
+                float[][] matrix =
+                [
+                    [1, 4, -1],
+                    [2, 5, 8],
+                    [1, 3, -3],
+                ];
+                float[] result = [4, 15, 1];
+                float[] expected = [1, 1, 1];
+                data.Add(matrix, result, expected);
+            }
+
+            {
+                float[][] matrix =
+                [
+                    [-1, 0, 0],
+                    [0, 1, 0],
+                    [0, 0, 1],
+                ];
+                float[] result = [1, 2, 3];
+                float[] expected = [-1, 2, 3];
+                data.Add(matrix, result, expected);
+            }
+
+            return data;
+        }
+    }
+}

tests/ImageSharp.Tests/Processing/Transforms/ProjectiveTransformTests.cs

@@ -55,6 +55,14 @@ public class ProjectiveTransformTests
         { TaperSide.Right, TaperCorner.RightOrBottom },
     };
 
+    public static readonly TheoryData<PointF, PointF, PointF, PointF> QuadDistortionData = new()
+    {
+        { new PointF(0, 0), new PointF(150, 0), new PointF(150, 150), new PointF(0, 150) }, // source == destination
+        { new PointF(25, 50), new PointF(210, 25), new PointF(140, 210), new PointF(15, 125) }, // Distortion
+        { new PointF(-50, -50), new PointF(200, -50), new PointF(200, 200), new PointF(-50, 200) }, // Scaling
+        { new PointF(150, 0), new PointF(150, 150), new PointF(0, 150), new PointF(0, 0) }, // Rotation
+    };
+
     public ProjectiveTransformTests(ITestOutputHelper output) => this.Output = output;
 
     [Theory]
@@ -93,6 +101,24 @@ public void Transform_WithTaperMatrix<TPixel>(TestImageProvider<TPixel> provider
         }
     }
 
+    [Theory]
+    [WithTestPatternImages(nameof(QuadDistortionData), 150, 150, PixelTypes.Rgba32)]
+    public void Transform_WithQuadDistortion<TPixel>(TestImageProvider<TPixel> provider, PointF topLeft, PointF topRight, PointF bottomRight, PointF bottomLeft)
+        where TPixel : unmanaged, IPixel<TPixel>
+    {
+        using (Image<TPixel> image = provider.GetImage())
+        {
+            ProjectiveTransformBuilder builder = new ProjectiveTransformBuilder()
+                .AppendQuadDistortion(topLeft, topRight, bottomRight, bottomLeft);
+
+            image.Mutate(i => i.Transform(builder));
+
+            FormattableString testOutputDetails = $"{topLeft}-{topRight}-{bottomRight}-{bottomLeft}";
+            image.DebugSave(provider, testOutputDetails);
+            image.CompareFirstFrameToReferenceOutput(TolerantComparer, provider, testOutputDetails);
+        }
+    }
+
     [Theory]
     [WithSolidFilledImages(100, 100, 0, 0, 255, PixelTypes.Rgba32)]
     public void RawTransformMatchesDocumentedExample<TPixel>(TestImageProvider<TPixel> provider)