diff --git a/src/ImageSharp/Common/Helpers/GaussianEliminationSolver.cs b/src/ImageSharp/Common/Helpers/GaussianEliminationSolver.cs new file mode 100644 index 0000000000..3bb9eb3956 --- /dev/null +++ b/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; + +/// +/// 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 +/// +/// The type of numbers used in the matrix and solution vector. Must implement the interface. +internal static class GaussianEliminationSolver + where TNumber : INumber +{ + /// + /// Solves the system of linear equations represented by the given matrix and result vector using Gaussian Elimination. + /// + /// The square matrix representing the coefficients of the linear equations. + /// The vector representing the constants on the right-hand side of the linear equations. + /// Thrown if the matrix is singular and cannot be solved. + /// + /// The matrix passed to this method must be a square matrix. + /// If the matrix is singular (i.e., has no unique solution), an will be thrown. + /// + 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]; + } + } + } +} diff --git a/src/ImageSharp/Common/Helpers/QuadDistortionHelper.cs b/src/ImageSharp/Common/Helpers/QuadDistortionHelper.cs new file mode 100644 index 0000000000..9fd6e84c75 --- /dev/null +++ b/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; + +/// +/// Provides helper methods for performing quad distortion transformations. +/// +internal static class QuadDistortionHelper +{ + /// + /// Computes the projection matrix for a quad distortion transformation. + /// + /// The source rectangle. + /// The top-left point of the distorted quad. + /// The top-right point of the distorted quad. + /// The bottom-right point of the distorted quad. + /// The bottom-left point of the distorted quad. + /// The computed projection matrix for the quad distortion. + /// + /// This method is based on the algorithm described in the following article: + /// https://blog.mbedded.ninja/mathematics/geometry/projective-transformations/ + /// + 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.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; + } +} diff --git a/src/ImageSharp/Processing/ProjectiveTransformBuilder.cs b/src/ImageSharp/Processing/ProjectiveTransformBuilder.cs index 0387adebb9..037a883291 100644 --- a/src/ImageSharp/Processing/ProjectiveTransformBuilder.cs +++ b/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 class ProjectiveTransformBuilder public ProjectiveTransformBuilder AppendTranslation(Vector2 position) => this.AppendMatrix(Matrix4x4.CreateTranslation(new Vector3(position, 0))); + /// + /// Prepends a quad distortion matrix using the specified corner points. + /// + /// The top-left corner point of the distorted quad. + /// The top-right corner point of the distorted quad. + /// The bottom-right corner point of the distorted quad. + /// The bottom-left corner point of the distorted quad. + /// The . + 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)); + + /// + /// Appends a quad distortion matrix using the specified corner points. + /// + /// The top-left corner point of the distorted quad. + /// The top-right corner point of the distorted quad. + /// The bottom-right corner point of the distorted quad. + /// The bottom-left corner point of the distorted quad. + /// The . + 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)); + /// /// Prepends a raw matrix. /// diff --git a/tests/ImageSharp.Tests/Common/GaussianEliminationSolverTest.cs b/tests/ImageSharp.Tests/Common/GaussianEliminationSolverTest.cs new file mode 100644 index 0000000000..85bc68b826 --- /dev/null +++ b/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.Solve(matrix, result); + + for (int i = 0; i < expected.Length; i++) + { + Assert.Equal(result[i], expected[i], 4); + } + } + + public static TheoryData MatrixTestData + { + get + { + TheoryData 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; + } + } +} diff --git a/tests/ImageSharp.Tests/Processing/Transforms/ProjectiveTransformTests.cs b/tests/ImageSharp.Tests/Processing/Transforms/ProjectiveTransformTests.cs index 4aec530364..293d39d663 100644 --- a/tests/ImageSharp.Tests/Processing/Transforms/ProjectiveTransformTests.cs +++ b/tests/ImageSharp.Tests/Processing/Transforms/ProjectiveTransformTests.cs @@ -55,6 +55,14 @@ public class ProjectiveTransformTests { TaperSide.Right, TaperCorner.RightOrBottom }, }; + public static readonly TheoryData 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 class ProjectiveTransformTests } } + [Theory] + [WithTestPatternImages(nameof(QuadDistortionData), 150, 150, PixelTypes.Rgba32)] + public void Transform_WithQuadDistortion(TestImageProvider provider, PointF topLeft, PointF topRight, PointF bottomRight, PointF bottomLeft) + where TPixel : unmanaged, IPixel + { + using (Image 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(TestImageProvider provider) diff --git a/tests/Images/External/ReferenceOutput/ProjectiveTransformTests/Transform_WithQuadDistortion_Rgba32_TestPattern150x150_PointF [ X=-50, Y=-50 ]-PointF [ X=200, Y=-50 ]-PointF [ X=200, Y=200 ]-PointF [ X=-50, Y=200 ].png b/tests/Images/External/ReferenceOutput/ProjectiveTransformTests/Transform_WithQuadDistortion_Rgba32_TestPattern150x150_PointF [ X=-50, Y=-50 ]-PointF [ X=200, Y=-50 ]-PointF [ X=200, Y=200 ]-PointF [ X=-50, Y=200 ].png new file mode 100644 index 0000000000..cac4e4ab24 --- /dev/null +++ b/tests/Images/External/ReferenceOutput/ProjectiveTransformTests/Transform_WithQuadDistortion_Rgba32_TestPattern150x150_PointF [ X=-50, Y=-50 ]-PointF [ X=200, Y=-50 ]-PointF [ X=200, Y=200 ]-PointF [ X=-50, Y=200 ].png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:fb2804f09a350d9889d57bec60972b6a60ac341f179926cbb8361e9809e47883 +size 33139 diff --git a/tests/Images/External/ReferenceOutput/ProjectiveTransformTests/Transform_WithQuadDistortion_Rgba32_TestPattern150x150_PointF [ X=0, Y=0 ]-PointF [ X=150, Y=0 ]-PointF [ X=150, Y=150 ]-PointF [ X=0, Y=150 ].png b/tests/Images/External/ReferenceOutput/ProjectiveTransformTests/Transform_WithQuadDistortion_Rgba32_TestPattern150x150_PointF [ X=0, Y=0 ]-PointF [ X=150, Y=0 ]-PointF [ X=150, Y=150 ]-PointF [ X=0, Y=150 ].png new file mode 100644 index 0000000000..2c4cb321b4 --- /dev/null +++ b/tests/Images/External/ReferenceOutput/ProjectiveTransformTests/Transform_WithQuadDistortion_Rgba32_TestPattern150x150_PointF [ X=0, Y=0 ]-PointF [ X=150, Y=0 ]-PointF [ X=150, Y=150 ]-PointF [ X=0, Y=150 ].png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0f468d2def55cda480d1adb077d8d087b4b15229d25bd39cd9e4fe7a203d6747 +size 1982 diff --git a/tests/Images/External/ReferenceOutput/ProjectiveTransformTests/Transform_WithQuadDistortion_Rgba32_TestPattern150x150_PointF [ X=150, Y=0 ]-PointF [ X=150, Y=150 ]-PointF [ X=0, Y=150 ]-PointF [ X=0, Y=0 ].png b/tests/Images/External/ReferenceOutput/ProjectiveTransformTests/Transform_WithQuadDistortion_Rgba32_TestPattern150x150_PointF [ X=150, Y=0 ]-PointF [ X=150, Y=150 ]-PointF [ X=0, Y=150 ]-PointF [ X=0, Y=0 ].png new file mode 100644 index 0000000000..cf48a0dcf8 --- /dev/null +++ b/tests/Images/External/ReferenceOutput/ProjectiveTransformTests/Transform_WithQuadDistortion_Rgba32_TestPattern150x150_PointF [ X=150, Y=0 ]-PointF [ X=150, Y=150 ]-PointF [ X=0, Y=150 ]-PointF [ X=0, Y=0 ].png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e5db33fb6cd26e7641370ece52710e3ec3ab2e23da1b4586616f2fa6ce0c4d46 +size 3030 diff --git a/tests/Images/External/ReferenceOutput/ProjectiveTransformTests/Transform_WithQuadDistortion_Rgba32_TestPattern150x150_PointF [ X=25, Y=50 ]-PointF [ X=210, Y=25 ]-PointF [ X=140, Y=210 ]-PointF [ X=15, Y=125 ].png b/tests/Images/External/ReferenceOutput/ProjectiveTransformTests/Transform_WithQuadDistortion_Rgba32_TestPattern150x150_PointF [ X=25, Y=50 ]-PointF [ X=210, Y=25 ]-PointF [ X=140, Y=210 ]-PointF [ X=15, Y=125 ].png new file mode 100644 index 0000000000..e0be7d88d5 --- /dev/null +++ b/tests/Images/External/ReferenceOutput/ProjectiveTransformTests/Transform_WithQuadDistortion_Rgba32_TestPattern150x150_PointF [ X=25, Y=50 ]-PointF [ X=210, Y=25 ]-PointF [ X=140, Y=210 ]-PointF [ X=15, Y=125 ].png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9a1ff90af73904b689c8ef1a63f0194c86ec46596a51c9c6cff7483d326b9968 +size 36739