From 45b2dee8bb164b45fd9bfb17acaa9c6fa6656b44 Mon Sep 17 00:00:00 2001 From: James Jackson-South Date: Mon, 12 Sep 2022 16:34:46 +1000 Subject: [PATCH 1/2] Add TryGetLinearlySeparableComponents and tests --- .../ConvolutionProcessorHelpers.cs | 104 +++++++++++++++++- .../ConvolutionProcessorHelpersTest.cs | 73 ++++++++++++ 2 files changed, 175 insertions(+), 2 deletions(-) create mode 100644 tests/ImageSharp.Tests/Processing/Processors/Convolution/ConvolutionProcessorHelpersTest.cs diff --git a/src/ImageSharp/Processing/Processors/Convolution/ConvolutionProcessorHelpers.cs b/src/ImageSharp/Processing/Processors/Convolution/ConvolutionProcessorHelpers.cs index ef25de6396..7dea3338ce 100644 --- a/src/ImageSharp/Processing/Processors/Convolution/ConvolutionProcessorHelpers.cs +++ b/src/ImageSharp/Processing/Processors/Convolution/ConvolutionProcessorHelpers.cs @@ -11,6 +11,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); @@ -18,9 +19,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; @@ -45,10 +48,12 @@ internal static float[] CreateGaussianBlurKernel(int size, float weight) /// /// Create a 1 dimensional Gaussian kernel using the Gaussian G(x) function /// + /// The kernel size. + /// The weight of the blur. /// The convolution kernel. internal static float[] CreateGaussianSharpenKernel(int size, float weight) { - var kernel = new float[size]; + float[] kernel = new float[size]; float sum = 0; @@ -85,5 +90,100 @@ 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..740bdddafd --- /dev/null +++ b/tests/ImageSharp.Tests/Processing/Processors/Convolution/ConvolutionProcessorHelpersTest.cs @@ -0,0 +1,73 @@ +// Copyright (c) Six Labors. +// Licensed under the Six Labors Split License. + +using System; +using SixLabors.ImageSharp.Processing.Processors.Convolution; +using Xunit; + +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; + } + } +} From 65088b362473c0d1544ba2b707cd38a7e3894140 Mon Sep 17 00:00:00 2001 From: James Jackson-South Date: Thu, 15 Sep 2022 22:12:29 +1000 Subject: [PATCH 2/2] Fix merge and namespace style --- .../ConvolutionProcessorHelpers.cs | 9 +- .../ConvolutionProcessorHelpersTest.cs | 97 +++++++++---------- 2 files changed, 54 insertions(+), 52 deletions(-) diff --git a/src/ImageSharp/Processing/Processors/Convolution/ConvolutionProcessorHelpers.cs b/src/ImageSharp/Processing/Processors/Convolution/ConvolutionProcessorHelpers.cs index ebda2eb079..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; diff --git a/tests/ImageSharp.Tests/Processing/Processors/Convolution/ConvolutionProcessorHelpersTest.cs b/tests/ImageSharp.Tests/Processing/Processors/Convolution/ConvolutionProcessorHelpersTest.cs index 740bdddafd..574c983710 100644 --- a/tests/ImageSharp.Tests/Processing/Processors/Convolution/ConvolutionProcessorHelpersTest.cs +++ b/tests/ImageSharp.Tests/Processing/Processors/Convolution/ConvolutionProcessorHelpersTest.cs @@ -1,73 +1,70 @@ // Copyright (c) Six Labors. // Licensed under the Six Labors Split License. -using System; using SixLabors.ImageSharp.Processing.Processors.Convolution; -using Xunit; -namespace SixLabors.ImageSharp.Tests.Processing.Processors.Convolution +namespace SixLabors.ImageSharp.Tests.Processing.Processors.Convolution; + +[GroupOutput("Convolution")] +public class ConvolutionProcessorHelpersTest { - [GroupOutput("Convolution")] - public class ConvolutionProcessorHelpersTest + [Theory] + [InlineData(3)] + [InlineData(5)] + [InlineData(9)] + [InlineData(22)] + [InlineData(33)] + [InlineData(80)] + public void VerifyGaussianKernelDecomposition(int radius) { - [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); + 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); + 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); + 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); + float[,] dotProduct = DotProduct(row, column); - for (int y = 0; y < column.Length; y++) + for (int y = 0; y < column.Length; y++) + { + for (int x = 0; x < row.Length; x++) { - for (int x = 0; x < row.Length; x++) - { - Assert.True(Math.Abs(matrix[y, x] - dotProduct[y, x]) < 0.0001F); - } + 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); + [Fact] + public void VerifyNonSeparableMatrix() + { + bool result = LaplacianKernels.LaplacianOfGaussianXY.TryGetLinearlySeparableComponents( + out float[] row, + out float[] column); - Assert.False(result); - Assert.Null(row); - Assert.Null(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]; + 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 x = 0; x < row.Length; x++) + { + for (int y = 0; y < column.Length; y++) { - for (int y = 0; y < column.Length; y++) - { - matrix[y, x] = row[x] * column[y]; - } + matrix[y, x] = row[x] * column[y]; } - - return matrix; } + + return matrix; } }