diff --git a/src/ImageSharp/Processing/Processors/Convolution/ConvolutionProcessorHelpers.cs b/src/ImageSharp/Processing/Processors/Convolution/ConvolutionProcessorHelpers.cs
index d2833fea9b..37adaf5402 100644
--- a/src/ImageSharp/Processing/Processors/Convolution/ConvolutionProcessorHelpers.cs
+++ b/src/ImageSharp/Processing/Processors/Convolution/ConvolutionProcessorHelpers.cs
@@ -9,6 +9,7 @@ internal static class ConvolutionProcessorHelpers
/// Kernel radius is calculated using the minimum viable value.
/// See .
///
+ /// The weight of the blur.
internal static int GetDefaultGaussianRadius(float sigma)
=> (int)MathF.Ceiling(sigma * 3);
@@ -16,9 +17,11 @@ internal static int GetDefaultGaussianRadius(float sigma)
/// Create a 1 dimensional Gaussian kernel using the Gaussian G(x) function.
///
/// The convolution kernel.
+ /// The kernel size.
+ /// The weight of the blur.
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
///
/// The convolution kernel.
+ /// The kernel size.
+ /// The weight of the blur.
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;
}
+
+ ///
+ /// 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, of size N and of size M.
+ /// This algorithm runs in O(NM).
+ ///
+ /// The input 2D matrix to analyze.
+ /// The resulting 1D row vector, if possible.
+ /// The resulting 1D column vector, if possible.
+ /// Whether or not was linearly separable.
+ public static bool TryGetLinearlySeparableComponents(this DenseMatrix 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;
+ }
}
diff --git a/tests/ImageSharp.Tests/Processing/Processors/Convolution/ConvolutionProcessorHelpersTest.cs b/tests/ImageSharp.Tests/Processing/Processors/Convolution/ConvolutionProcessorHelpersTest.cs
new file mode 100644
index 0000000000..574c983710
--- /dev/null
+++ b/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 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 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;
+ }
+}