Add TryGetLinearlySeparableComponents and tests

src/ImageSharp/Processing/Processors/Convolution/ConvolutionProcessorHelpers.cs

@@ -9,16 +9,19 @@ internal static class ConvolutionProcessorHelpers
     /// Kernel radius is calculated using the minimum viable value.
     /// See <see href="http://chemaguerra.com/gaussian-filter-radius/"/>.
     /// </summary>
+    /// <param name="sigma">The weight of the blur.</param>
     internal static int GetDefaultGaussianRadius(float sigma)
         => (int)MathF.Ceiling(sigma * 3);
 
     /// <summary>
     /// Create a 1 dimensional Gaussian kernel using the Gaussian G(x) function.
     /// </summary>
     /// <returns>The convolution kernel.</returns>
+    /// <param name="size">The kernel size.</param>
+    /// <param name="weight">The weight of the blur.</param>
     internal static float[] CreateGaussianBlurKernel(int size, float weight)
     {
-        var kernel = new float[size];
+        float[] kernel = new float[size];
 
         float sum = 0F;
         float midpoint = (size - 1) / 2F;
@@ -44,9 +47,11 @@ internal static float[] CreateGaussianBlurKernel(int size, float weight)
     /// Create a 1 dimensional Gaussian kernel using the Gaussian G(x) function
     /// </summary>
     /// <returns>The convolution kernel.</returns>
+    /// <param name="size">The kernel size.</param>
+    /// <param name="weight">The weight of the blur.</param>
     internal static float[] CreateGaussianSharpenKernel(int size, float weight)
     {
-        var kernel = new float[size];
+        float[] kernel = new float[size];
 
         float sum = 0;
 
@@ -83,4 +88,99 @@ internal static float[] CreateGaussianSharpenKernel(int size, float weight)
 
         return kernel;
     }
+
+    /// <summary>
+    /// Checks whether or not a given NxM matrix is linearly separable, and if so, it extracts the separable components.
+    /// These would be two 1D vectors, <paramref name="row"/> of size N and <paramref name="column"/> of size M.
+    /// This algorithm runs in O(NM).
+    /// </summary>
+    /// <param name="matrix">The input 2D matrix to analyze.</param>
+    /// <param name="row">The resulting 1D row vector, if possible.</param>
+    /// <param name="column">The resulting 1D column vector, if possible.</param>
+    /// <returns>Whether or not <paramref name="matrix"/> was linearly separable.</returns>
+    public static bool TryGetLinearlySeparableComponents(this DenseMatrix<float> matrix, out float[] row, out float[] column)
+    {
+        int height = matrix.Rows;
+        int width = matrix.Columns;
+
+        float[] tempX = new float[width];
+        float[] tempY = new float[height];
+
+        // This algorithm checks whether the input matrix is linearly separable and extracts two
+        // 1D components if possible. Note that for a given NxM matrix that is linearly separable,
+        // there exists an infinite number of possible solutions to the system of linear equations
+        // representing the possible 1D components that can produce the input matrix as a product.
+        // Let's assume we have a 3x3 input matrix to describe the logic. We have the following:
+        //
+        //     | m11, m12, m13 |                                                  | c1 |
+        // M = | m21, m22, m23 |, and we want to find: R = | r1, r2, r3 | and C = | c2 |.
+        //     | m31, m32, m33 |                                                  | c3 |
+        //
+        // We essentially get the following system of linear equations to solve:
+        //
+        //   / a11 = r1c1
+        //   | a12 = r2c1
+        //   | a13 = r3c1
+        //   | a21 = r1c2                   a11     a12     a13       a11     a12     a13
+        //   / a22 = r2c2, which gives us: ----- = ----- = ----- and ----- = ----- = -----.
+        //   \ a23 = r3c2                   a21     a22     a23       a31     a32     a33
+        //   | a31 = r1c3
+        //   | a32 = r2c3
+        //   \ a33 = r3c3
+        //
+        // As we said, there are infinite solutions to this problem (provided the input matrix is in
+        // fact linearly separable), but we can look at the equalities above to find a way to define
+        // one specific solution that is very easy to calculate (and that is equivalent to all others
+        // anyway). In particular, we can see that in order for it to be linearly separable, the matrix
+        // needs to have each row linearly dependent on each other. That is, its rank is just 1. This
+        // means that we can express the whole matrix as a function of one row vector (any of the rows
+        // in the matrix), and a column vector that represents the ratio of each element in a given column
+        // j with the corresponding j-th item in the reference row. This same procedure extends naturally
+        // to lineary separable 2D matrices of any size, too. So we end up with the following generalized
+        // solution for a matrix M of size NxN (or MxN, that works too) and the R and C vectors:
+        //
+        //     | m11, m12, m13, ..., m1N |                                       | m11/m11 |
+        //     | m21, m22, m23, ..., m2N |                                       | m21/m11 |
+        // M = | m31, m32, m33, ..., m3N |, R = | m11, m12, m13, ..., m1N |, C = | m31/m11 |.
+        //     | ...  ...  ...  ...  ... |                                       |   ...   |
+        //     | mN1, mN2, mN3, ..., mNN |                                       | mN1/m11 |
+        //
+        // So what this algorithm does is just the following:
+        //   1) It calculates the C[i] value for each i-th row.
+        //   2) It checks that every j-th item in the row respects C[i] = M[i, j] / M[0, j]. If this is
+        //      not true for any j-th item in any i-th row, then the matrix is not linearly separable.
+        //   3) It sets items in R and C to the values detailed above if the validation passed.
+        for (int y = 1; y < height; y++)
+        {
+            float ratio = matrix[y, 0] / matrix[0, 0];
+
+            for (int x = 1; x < width; x++)
+            {
+                if (Math.Abs(ratio - (matrix[y, x] / matrix[0, x])) > 0.0001f)
+                {
+                    row = null;
+                    column = null;
+
+                    return false;
+                }
+            }
+
+            tempY[y] = ratio;
+        }
+
+        // The first row is used as a reference, to the ratio is just 1
+        tempY[0] = 1;
+
+        // The row component is simply the reference row in the input matrix.
+        // In this case, we're just using the first one for simplicity.
+        for (int x = 0; x < width; x++)
+        {
+            tempX[x] = matrix[0, x];
+        }
+
+        row = tempX;
+        column = tempY;
+
+        return true;
+    }
 }

tests/ImageSharp.Tests/Processing/Processors/Convolution/ConvolutionProcessorHelpersTest.cs

@@ -0,0 +1,70 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Processing.Processors.Convolution;
+
+namespace SixLabors.ImageSharp.Tests.Processing.Processors.Convolution;
+
+[GroupOutput("Convolution")]
+public class ConvolutionProcessorHelpersTest
+{
+    [Theory]
+    [InlineData(3)]
+    [InlineData(5)]
+    [InlineData(9)]
+    [InlineData(22)]
+    [InlineData(33)]
+    [InlineData(80)]
+    public void VerifyGaussianKernelDecomposition(int radius)
+    {
+        int kernelSize = (radius * 2) + 1;
+        float sigma = radius / 3F;
+        float[] kernel = ConvolutionProcessorHelpers.CreateGaussianBlurKernel(kernelSize, sigma);
+        DenseMatrix<float> matrix = DotProduct(kernel, kernel);
+
+        bool result = matrix.TryGetLinearlySeparableComponents(out float[] row, out float[] column);
+
+        Assert.True(result);
+        Assert.NotNull(row);
+        Assert.NotNull(column);
+        Assert.Equal(row.Length, matrix.Rows);
+        Assert.Equal(column.Length, matrix.Columns);
+
+        float[,] dotProduct = DotProduct(row, column);
+
+        for (int y = 0; y < column.Length; y++)
+        {
+            for (int x = 0; x < row.Length; x++)
+            {
+                Assert.True(Math.Abs(matrix[y, x] - dotProduct[y, x]) < 0.0001F);
+            }
+        }
+    }
+
+    [Fact]
+    public void VerifyNonSeparableMatrix()
+    {
+        bool result = LaplacianKernels.LaplacianOfGaussianXY.TryGetLinearlySeparableComponents(
+            out float[] row,
+            out float[] column);
+
+        Assert.False(result);
+        Assert.Null(row);
+        Assert.Null(column);
+    }
+
+    private static DenseMatrix<float> DotProduct(float[] row, float[] column)
+    {
+        float[,] matrix = new float[column.Length, row.Length];
+
+        for (int x = 0; x < row.Length; x++)
+        {
+            for (int y = 0; y < column.Length; y++)
+            {
+                matrix[y, x] = row[x] * column[y];
+            }
+        }
+
+        return matrix;
+    }
+}